1 Problématique

Dans le secteur de l’assurance, en particulier dans les produits d’assurance automobile, la fidélité des clients est un enjeu majeur pour la pérennité des entreprises. La concurrence accrue dans ce secteur, la numérisation des services, ainsi que l’évolution des besoins des consommateurs ont poussé les compagnies d’assurance à réévaluer régulièrement leur capacité à retenir les clients. Ce phénomène de départ des assurés, également appelé churn ou désabonnement, est particulièrement critique puisqu’il implique non seulement une perte de revenu direct, mais aussi un impact sur la réputation et les marges de la compagnie.

Alors on pourrait se demander quels sont les facteurs socio-économiques, géographiques et comportementaux qui influencent la durée de fidélité des clients avant qu’ils ne résilient leur contrat d’assurance automobile ?

2 Présentation des données

L’objectif de cette étude est d’analyser le taux de désabonnement des assurés d’une compagnie d’assurance automobile afin de mieux comprendre les facteurs influençant le comportement de résiliation des contrats. En effet, une analyse approfondie de la durée de vie des clients, c’est-à-dire du temps écoulé entre la souscription à un contrat d’assurance et la résiliation, permet de mettre en lumière des tendances et des phénomènes récurrents, tels que la vulnérabilité des clients à la concurrence après une certaine période. L’analyse se concentrera sur les facteurs suivants :

Caractéristiques démographiques : Âge, état civil, présence d’enfants, niveau d’éducation (diplôme universitaire), etc. Facteurs économiques : Revenu, valeur marchande du logement, statut de propriétaire. Comportements liés à l’assurance : Ancienneté (tenure) dans le contrat, montant de la prime annuelle, qualité du crédit. Facteurs géographiques : Région, latitude, longitude, ville et état de résidence.

Nous avons effectuer des modifications au jeu de données initial qui vise à réduire le déséquilibre de classes dans la variable cible Churn et à se concentrer sur une partie plus intéressante des clients c’est-à-dire la partie où on observe plus de désabonnements et à des dates beaucoup plus variées. En outre, des ajustements ont été faits pour faciliter l’analyse des clients activés en 2020 (date d’effet du contrat) avec un seuil d’ancienneté significatif.

Pour équilibrer les classes de la variable cible, 90 % des observations où la variable Churn est égale à 0 ont été supprimées. Cela a été réalisé en identifiant les indices des lignes concernées, en sélectionnant aléatoirement les lignes à supprimer, puis en retirant ces lignes du jeu de données. Cette étape vise à atténuer le déséquilibre de la distribution des classes dans les données.

La colonne cust_orig_date, représentant la date d’origine des clients, a été convertie au format Date. À partir de cette colonne, une nouvelle variable start_year a été extraite, correspondant à l’année d’activation du client. Cette transformation permet de réaliser des analyses basées sur le temps, en particulier pour filtrer des cohortes temporelles spécifiques.

Le jeu de données a été restreint aux clients activés en 2020. Cela a permis de concentrer l’analyse sur un sous-groupe homogène et temporellement pertinent, facilitant une comparaison cohérente entre les individus de cette cohorte.

La variable start_year a été convertie en format numérique afin de permettre des manipulations mathématiques ou statistiques si nécessaire. La variable Churn a été transformée en facteur, ce qui est essentiel pour effectuer des analyses catégoriques ou pour entraîner des modèles adaptés à des variables qualitatives.

Pour se concentrer sur des clients fidèles ou à risque élevé, seuls ceux ayant une ancienneté (days_tenure) supérieure ou égale à 700 jours ont été conservés. Ce filtrage permet de recentrer l’analyse sur les clients ayant une contribution significative ou nécessitant une intervention particulière.

Une fois les données filtrées pour les clients ayant days_tenure >= 700, la valeur de cette variable a été recalculée en soustrayant 700. Cela permet de repositionner l’ancienneté par rapport à ce seuil et de simplifier l’interprétation des analyses ultérieures.

2.1 Description des variables

  • Variables qualitatives :
Variables Description
individual_id Identifiant unique de l’individu
address_id Identifiant de l’adresse associée à l’individu
cust_orig_date Date d’activation de l’abonnement client
date_of_birth Date de naissance de l’individu
city Ville de résidence
state État de résidence
county Comté de résidence
has_children Indicateur binaire : 1 si l’individu a des enfants, 0 sinon
marital_status Statut marital de l’individu (célibataire, marié, etc.)
home_owner Propriétaire du logement : 1 si propriétaire, 0 sinon
college_degree Diplôme universitaire : 1 si l’individu a un diplôme, 0 sinon
good_credit Indicateur de bon crédit : 1 si l’individu a un bon crédit, 0 sinon
acct_suspd_date Date de suspension du compte
Churn 1 si l’individu a résilié son abonnement, 0 sinon
  • Variables quantitatives :
Variables Description
curr_ann_amt Montant annuel actuel payé par le client
days_tenure Nombre de jours depuis que le client a souscrit à l’abonnement
age_in_years Âge de l’individu en années
latitude Latitude de la résidence de l’individu
longitude Longitude de la résidence de l’individu
income Revenu annuel de l’individu
length_of_residence Durée de résidence de l’individu à l’adresse actuelle
home_market_value Valeur marchande estimée du logement de l’individu
  • Statistiques descriptives
##  individual_id         address_id         curr_ann_amt      days_tenure   
##  Min.   :2.213e+11   Min.   :5.213e+11   Min.   :  32.99   Min.   :  1.0  
##  1st Qu.:2.213e+11   1st Qu.:5.213e+11   1st Qu.: 775.56   1st Qu.: 84.0  
##  Median :2.213e+11   Median :5.213e+11   Median : 932.06   Median :160.0  
##  Mean   :2.213e+11   Mean   :5.213e+11   Mean   : 934.32   Mean   :170.4  
##  3rd Qu.:2.213e+11   3rd Qu.:5.213e+11   3rd Qu.:1093.94   3rd Qu.:266.0  
##  Max.   :2.213e+11   Max.   :5.213e+11   Max.   :1815.52   Max.   :366.0  
##                                                                           
##  cust_orig_date        age_in_years   date_of_birth         latitude    
##  Min.   :2020-01-01   Min.   :24.00   Length:11068       Min.   :32.17  
##  1st Qu.:2020-04-10   1st Qu.:37.00   Class :character   1st Qu.:32.71  
##  Median :2020-07-25   Median :43.00   Mode  :character   Median :32.84  
##  Mean   :2020-07-14   Mean   :43.46                      Mean   :32.85  
##  3rd Qu.:2020-10-09   3rd Qu.:49.00                      3rd Qu.:32.99  
##  Max.   :2020-12-31   Max.   :55.00                      Max.   :33.51  
##                                                          NA's   :2114   
##    longitude          city              state              county         
##  Min.   :-97.96   Length:11068       Length:11068       Length:11068      
##  1st Qu.:-97.16   Class :character   Class :character   Class :character  
##  Median :-96.90   Mode  :character   Mode  :character   Mode  :character  
##  Mean   :-96.94                                                           
##  3rd Qu.:-96.71                                                           
##  Max.   :-96.12                                                           
##  NA's   :2114                                                             
##      income        has_children    length_of_residence marital_status    
##  Min.   :  5000   Min.   :0.0000   Min.   : 0.000      Length:11068      
##  1st Qu.: 42500   1st Qu.:0.0000   1st Qu.: 2.000      Class :character  
##  Median : 70000   Median :1.0000   Median : 5.000      Mode  :character  
##  Mean   : 74029   Mean   :0.6183   Mean   : 5.704                        
##  3rd Qu.: 87500   3rd Qu.:1.0000   3rd Qu.: 8.000                        
##  Max.   :250000   Max.   :1.0000   Max.   :15.000                        
##                                                                          
##  home_market_value    home_owner     college_degree    good_credit    
##  Length:11068       Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
##  Class :character   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:1.0000  
##  Mode  :character   Median :1.0000   Median :0.0000   Median :1.0000  
##                     Mean   :0.7481   Mean   :0.2902   Mean   :0.8154  
##                     3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1.0000  
##                     Max.   :1.0000   Max.   :1.0000   Max.   :1.0000  
##                                                                       
##  acct_suspd_date    Churn      start_year  
##  Length:11068       0:6209   Min.   :2020  
##  Class :character   1:4859   1st Qu.:2020  
##  Mode  :character            Median :2020  
##                              Mean   :2020  
##                              3rd Qu.:2020  
##                              Max.   :2020  
##

2.2 Analyse exploratoire

3 Estimation non paramétrique du modèle sans covariable en présence ou non de censure

L’analyse de survie, aussi appelée fiabilité, est utilisée pour étudier et modéliser le temps écoulé T jusqu’à un événement spécifique, tel que le défaut de paiement d’un emprunteur ou la résiliation d’un contrat d’assurance. Cette approche permet de prendre en compte des données censurées, c’est-à-dire des cas où l’événement n’a pas encore eu lieu lorsque l’étude se termine.

Dans notre cas,

T : Temps écoulé avant la résilliation d’un contrat qui correspond à la durée de vie du contrat d’assurance avant qu’il ne soit résilié (désabonnement).

S(t) : La fonction de survie représente la probabilité qu’un client n’ait pas résilié son contrat avant un moment donné t. Elle indique donc la proportion d’assurés qui sont encore sous contrat à un temps t. Une valeur de S(t) proche de 1 signifie que peu de résiliations ont eu lieu, tandis qu’une valeur proche de 0 indique que presque tous les assurés ont résilié. \[S(t) = P(T > t) = 1 - F(t)\]

h(t) : La fonction de risque instantané h(t), correspond au risque instantané de résiliation d’un contrat au temps t, sachant que le client est encore sous contrat à ce moment-là. Elle est définie comme la probabilité conditionnelle de résiliation dans un intervalle de temps très court autour de t, donné que le client a survécu jusqu’à t. Elle nous renseigne sur la vitesse à laquelle les résiliations surviennent à un instant donné. \[h(t) = \frac{f(t)}{S(t)}\]

H(t) : La fonction de risque cumulé H(t) mesure l’accumulation du risque de résiliation jusqu’à un moment t. Elle représente la somme des risques instantanés sur la période de temps écoulée. Cette fonction est utile pour quantifier l’intensité totale du risque jusqu’à un moment donné, offrant une vue d’ensemble du risque global auquel les assurés ont été exposés pendant la durée étudiée. \[H(t) = \int_{0}^{t} h(x) \, \mathrm{d}x = - ln(S(t))\]

Nous avons une colonne Churn dans notre table de données qui représente notre variable censure. Cette variable prend la valeur 0 si l’évènement d’intérêt (dans notre cas la résilliation du contrat) n’a pas eu lieu pendant la période d’observation et 1 s’il a eu lieu.

3.1 Fonction de survie S(t)

Sans Censure

En absence de censure, l’estimateur de la fonction de survie est donnée par : \[\hat{S}(t) = \frac{1}{n}\sum_{i=1}^n 1_{t_{i}>t}\]

## Call: survfit(formula = surv ~ 1, data = data_1, conf.type = "log-log")
## 
##          n events median 0.95LCL 0.95UCL
## [1,] 11068   4859    287     282     291

Avec censure

L’estimateur de notre fonction de survie S est definit par l’estimateur de Kaplan Meier : \[\hat{S_{K M}}(t_i) = S(t_{i-1})(1-\frac{d_i}{n_i}) = \prod_{t_{i}< t}(1 - \frac{d_{i}}{n_{i}}) \]

avec :

  • \(\hat{S_{K M}}(t_i)\) : la probabilité qu’il n’y ait pas de résilliation à l’instant t
  • \(t_i\) : Les temps où des événements (résiliations) se produisent.
  • \(d_i\) : Le nombre d’événements (résiliations) observés à \(t_i\).
  • \(n_i\) : Le nombre de clients “en risque” juste avant \(t_i\) (c’est-à-dire ceux qui n’ont pas encore résilié ou ne sont pas censurés).

Démonstration de l’estimateur de Kaplan-Meier

L’objectif est de démontrer comment l’estimateur de Kaplan-Meier est obtenu pour la fonction de survie.

  • Étape 1 : Définition récursive de l’estimateur de Kaplan-Meier

L’estimateur de Kaplan-Meier repose sur une mise à jour de la fonction de survie à chaque instant \(t_i\), où un événement (résiliation) survient. La probabilité de survie entre deux instants \(t_{i-1}\) et \(t_i\) est donnée par la probabilité de survie jusqu’à \(t_{i-1}\), multipliée par la probabilité conditionnelle de ne pas résilier à l’instant \(t_i\).

Ainsi, à chaque instant \(t_i\), la fonction de survie est mise à jour selon la formule :

\[ S(t_i) = S(t_{i-1}) \times \left( 1 - \frac{d_i}{n_i} \right) \]

  • Étape 2 : Expansion de la formule récursive

La relation est itérative, ce qui signifie que pour chaque instant \(t_i\), la fonction de survie est le produit des probabilités conditionnelles jusqu’à cet instant. Prenons quelques exemples pour comprendre comment cela fonctionne :

Pour le premier événement \(t_1\) :

\[ S(t_1) = 1 - \frac{d_1}{n_1} \]

Pour le second événement \(t_2\) :

\[ S(t_2) = S(t_1) \times \left( 1 - \frac{d_2}{n_2} \right) \] \[ S(t_2) = \left( 1 - \frac{d_1}{n_1} \right) \times \left( 1 - \frac{d_2}{n_2} \right) \]

Et ainsi de suite pour \(t_3\) :

\[ S(t_3) = S(t_2) \times \left( 1 - \frac{d_3}{n_3} \right) \] \[ S(t_3) = \left( 1 - \frac{d_1}{n_1} \right) \times \left( 1 - \frac{d_2}{n_2} \right) \times \left( 1 - \frac{d_3}{n_3} \right) \]

  • Étape 3 : Formulation générale

En généralisant cette relation à tous les instants \(t_i\), on peut écrire la fonction de survie sous forme de produit pour tout \(t\) comme suit :

\[ \hat{S_{KM}}(t) = \prod_{t_i \leq t} \left( 1 - \frac{d_i}{n_i} \right) \]

Cette formule signifie que la probabilité de survie à un temps \(t\) est le produit des probabilités de survie conditionnelles pour tous les instants où des événements se sont produits avant \(t\). Cela tient compte des résiliations successives et du nombre d’individus à risque à chaque instant.

Ainsi, la fonction de survie estimée par Kaplan-Meier repose sur cette accumulation des probabilités de ne pas subir l’événement (résiliation) à chaque instant, jusqu’à \(t\).

## Call: survfit(formula = new_surv ~ 1, data = data_1, conf.type = "log-log")
## 
##          n events median 0.95LCL 0.95UCL
## [1,] 11068   6209    252     248     257

Les propriétés de l’estimateur de Kaplan-Meier

Biais et convergence

L’estimateur de Kaplan-Meier est un estimateur légèrement biaisé: \[E(\hat{S_{K M}}(t)) < {S}(t)\] où \(E\) désigne l’espérance de l’estimateur de Kaplan-Meier.

Il est en revanche convergent : \[\forall\epsilon > 0, \; \lim_{n\to\infty} \mathbf{P}\left(|\hat{S_{K M}}(t)-{S}(t)|\ge\epsilon\right) = 0\]

Il est donc asymptotiquement sans biais : \[\lim_{n\to\infty} E(\hat{S_{K M}}(t)) = {S}(t)\]

Variance

  • la variance de l’estimateur est donnée par : \[\hat{Var}(\hat{S_{K M}(t)}) = \hat{S_{K M}}(t)²\sum_{t_i<t}\frac{d_i}{n_i(n_i-d_i)}\]

Normalité asymptotique

  • Propriété : L’estimateur de Kaplan-Meier \(\hat{S_{K M}}\) suit asymptotiquement une loi Normale, cette propriété permet de déterminer un intervalle de confiance à partir de la variance de cet estimateur.

En tout point de continuité de S :

\[\sqrt{n}(\hat{S_{KM}}(t) - S(t))\xrightarrow{L} N(0,I(S(t))^{-1})\]

On peut donc faire une approximation par la loi normale et ainsi obtenir un intervalle de confiance à 5% :

\[IC(\alpha) = \left[\hat{S_{K M}}(t) \pm z_{\frac{\alpha}{2}}\sqrt{\hat{Var}(\hat{S_{K M}(t)})} \right]\]

On détermine la variance de l’estimateur de Kaplan-Meyer pour la fonction de survie avec l’estimateur de Greenwood qui permet d’obtenir l’intervalle de confiance.

Nom des paramètres Formule Description
Estimateur de Greenwood \(σ_{KM}^2 = var(\hat{S}_{KM})\) C’est la variance de l’estimateur de Kaplan-Meyer pour la fonction de survie
  • Sans Censure
## Limite inférieure sans censure: 0.9990242 0.9980717 0.9973989 0.9971788 0.9951533 0.9943213
## Limite supérieure sans censure : 0.9998915 0.9993969 0.9989831 0.9988406 0.9974218 0.9968014
  • Avec Censure
## Limite inférieure AVEC censure: 0.9983027 0.9968559 0.9957894 0.9955607 0.9938822 0.9926173
## Limite supérieure AVEC censure : 0.9995289 0.9986253 0.9978825 0.9981112 0.9967117 0.9958

Représentation graphique de la fonction de survie

  • Sans censure : On suit le taux de survie des clients sans prendre en compte les données censurées (c’est-à-dire les contrats qui n’ont pas encore été résiliés au moment de la collecte des données). La courbe de survie présente deux chutes significatives : une première autour de 1 mois an d’ancienneté et une autre autour d’un peu plus de 3 mois. Ces baisses indiquent que des résiliations importantes ont eu lieu à ces moments, réduisant la probabilité de survie de façon marquée.

  • Avec censure : On tient compte ici des contrats qui sont toujours actifs à la fin de la période d’étude. On observe des baisses de la fonction de survie similaires à celles de la courbe sans censure. L’inclusion de la censure rend l’analyse plus réaliste, car elle montre que certains clients sont encore sous contrat au-delà d’une année, ce qui pourrait être interprété comme un signe de fidélité.

3.2 Fonction de risque instantané h(t)

Représentation graphique de la fonction de risque

Nous constatons que la fonction risque instantané SANS censure et AVEC censure sont quasiment les mêmes.

La similarité des deux graphes de la fonction de risque indique que la censure n’affecte pas significativement la dynamique du risque de désabonnement dans notre ensemble de données.

3.3 Fonction de risque cumulé H(t)

Lorsqu’il y a présence de données censurées, 2 estimateurs de H(t) sont proposés :

3.3.1 Estimation par Breslow

SANS censure AVEC censure
\(\hat{H}_{Breslow}(t) = -log~(\hat{S}(t))\) \(\hat{H}_{Breslow}(t) = -log~(\hat{S_{K M}}(t)) = -\sum_{t_i \leq t} log \left(1 - \frac{d_i}{n_i}\right)\)

La variance de l’estimateur de Breslow est donnée par : \[\hat{Var}(\hat{H}_{Breslow}(t)) = \sum_{t_i \le t} \frac{d_i}{n_i(n_i - d_i)}\]

  • Représentation graphique

On observe que le risque cumulé sans censure augmente plus rapidement qu’avec censure, en utilisant l’estimateur de Breslow.

3.3.2 Estimation par Nelson-Aalen

Il est basé sur la relation : \(\int_0^t H(x) dx\) et est definit par: \[\hat{H}_{NA}(t) = \sum_{i=1}^n\frac{d_i}{n_i} \times \mathbb{I_{(t_i <= t)}}\]

La variance de l’estimateur de Nelson-Aalen est : \[\hat{Var}(\hat{H}_{NA}(t)) = \sum_{t_i \leq t} \frac{d_i}{n_i²}\]

  • Représentation graphqiue

  • En présence et en absence de censure, on observe un léger écart en fin de période avec l’estimation par Nelson-Aalen.

  • Dans tous les cas le risque cumulé augmente progressivement. Ce qui veut dire que le risque de résilliation de contrat devient plus important à chaque fois qu’on avance au cours du temps qui pourrait être expliqué par d’autres facteurs représentés par certains variables dans notre jeu de données.

3.3.3 Comparaison : estimateurs de Nelson-Aalen et Breslow

Les deux estimateurs sont quasi-similaires pour estimer la fonction de risque cumulé H(t).

4 Estimation et test paramétrique du modèle sans covariable en présence ou non de censure

Hypothèses :

  • H0 : Les durées d’abonnement (temps avant résiliation) suivent une distribution théorique donnée (exponentielle, Weibull, gamma et log-normal).

  • H1 : Les durées d’abonnement (temps avant résiliation) ne suivent pas la distribution théorique donnée.

Tableau des modèles :

Modèle Formule de répartition \(F(t)\) Fonction de survie \(S(t)\)
Weibull \(F(t) = 1 - e^{-(\frac{t - \gamma}{\mu})^k}, \quad \text{si } t \geqslant \gamma\) \(S(t) = e^{-(\frac{t - \gamma}{\mu})^k}, \quad \text{si } t \geqslant \gamma\)
Exponentielle \(F(t) = 1 - e^{-\lambda t}\) \(S(t) = e^{-\lambda t}\)
Gamma \(F(t) = \frac{\int_0^t x^{\alpha - 1} e^{-x} \, dx}{\Gamma(\alpha)}\) \(S(t) = 1 - \frac{\int_0^t x^{\alpha - 1} e^{-x} \, dx}{\Gamma(\alpha)}\)
Log-Normale \(F(t) = \Phi \left( \frac{\ln(t) - \mu}{\sigma} \right)\) \(S(t) = 1 - \Phi \left( \frac{\ln(t) - \mu}{\sigma} \right)\)
  • \(\Phi\) : Fonction de répartition cumulative de la loi normale standard (fonction “normale cumulative”).
  • \(\mu\) et \(\sigma\) : Paramètres de la loi log-normale, représentant respectivement la moyenne et l’écart-type de la variable \(\ln(T)\) (le logarithme du temps).
  • \(\lambda\) : Taux de résiliation pour la loi exponentielle.
  • \(\gamma\), \(\mu\), et \(k\) : Paramètres de la loi de Weibull. \(\gamma\) est le décalage, \(\mu\) est l’échelle, et \(k\) est la forme de la distribution.
  • \(\alpha\) : Paramètre de forme pour la loi gamma.
  • \(\Gamma(\alpha)\) : Fonction gamma d’Euler, utilisée dans le calcul de la loi gamma, définie comme une généralisation de la factorielle pour les réels et les complexes.

Le tracé d’un objet de classe fitdist fournit quatre tracés classiques de qualité d’ajustement présentés sur la figure ci-dessous:

un diagramme de densité représentant la fonction de densité de la distribution ajustée avec l’histogramme de la distribution empirique,

un graphique CDF de la distribution empirique et de la distribution ajustée,

un graphique Q-Q représentant les quantiles empiriques (axe y) par rapport aux quantiles théoriques (axe x)

un graphique P-P représentant la fonction de distribution empirique évaluée à chaque point de données (axe y) par rapport à la fonction de distribution ajustée (axe x)

4.1 Loi weibull

La fonction de répartition de la loi de Weibull est :
\(F(t) = 0, si : t < \gamma\)
avec \(\beta\), \(\mu > 0\), \(\gamma \geqslant 0\).

où :

  • \(\gamma\) : Décalage à l’origine du temps
  • \(\beta\) : Paramètre de forme (donne le type de défaillance étudié)
  • \(\mu\) : Paramètre d’échelle

Weibull sans censure

## Fitting of the distribution ' weibull ' by maximum likelihood 
## Parameters:
##         estimate Std. Error
## shape   1.470181 0.01174164
## scale 187.027663 1.26587759

## AIC Weibull: 133876.4
## BIC Weibull: 133891

Weibull avec censure

## Call:
## survreg(formula = Surv(days_tenure, censor_status) ~ 1, data = data_1, 
##     dist = "weibull")
## 
## Coefficients:
## (Intercept) 
##    5.796473 
## 
## Scale= 0.6920544 
## 
## Loglik(model)= -33423.3   Loglik(intercept only)= -33423.3
## n= 11068

Test weibull

  • Sans censure
## 
## ### Résultats du modèle Weibull ###
## 
## Test de Kolmogorov-Smirnov :
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  positive_data
## D = 0.31168, p-value < 2.2e-16
## alternative hypothesis: two-sided
## 
## Test de Cramer-von Mises :
## 
##  Cramer-von Mises test of goodness-of-fit
##  Null hypothesis: Weibull distribution
##  with parameters shape = 1.44497315484171, scale = 329.136531910119
##  Parameters assumed to be fixed
## 
## data:  positive_data
## omega2 = 434.54, p-value < 2.2e-16
## 
## Test d'Anderson-Darling :
## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: Weibull distribution
##  with parameters shape = 1.44497315484171, scale = 329.136531910119
##  Parameters assumed to be fixed
## 
## data:  positive_data
## An = 2407.9, p-value = 5.421e-08
  • Avec censure
##              est       L95%       U95%         se
## shape   1.444973   1.411177   1.479579 0.01744828
## scale 329.136537 322.335681 336.080882 3.50623847
## 
## ### Résultats du modèle Weibull ###
## 
## Test de Kolmogorov-Smirnov :
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  positive_data
## D = 0.31168, p-value < 2.2e-16
## alternative hypothesis: two-sided
## 
## Test de Cramer-von Mises :
## 
##  Cramer-von Mises test of goodness-of-fit
##  Null hypothesis: Weibull distribution
##  with parameters shape = 1.44497314671028, scale = 329.136537251605
##  Parameters assumed to be fixed
## 
## data:  positive_data
## omega2 = 434.54, p-value < 2.2e-16
## 
## Test d'Anderson-Darling :
## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: Weibull distribution
##  with parameters shape = 1.44497314671028, scale = 329.136537251605
##  Parameters assumed to be fixed
## 
## data:  positive_data
## An = 2407.9, p-value = 5.421e-08

4.2 Loi exponentielle

La loi exponentielle est très utilisée dans les études d’analyse de survie, puisqu’elle est “sans mémoire”.

Exponentielle sans censure

## Fitting of the distribution ' exp ' by matching moments 
## Parameters : 
##         estimate  Std. Error
## rate 0.005868822 0.005868822
## Loglikelihood:  -67936.51   AIC:  135875   BIC:  135882.3

##        rate    
##   5.868822e-03 
##  (5.578486e-05)
## AIC exponentielle: 135875
## BIC exponentielle: 135882.3

Exponentielle avec censure

## Call:
## survreg(formula = Surv(days_tenure, censor_status) ~ 1, data = data_1, 
##     dist = "exponential")
## 
## Coefficients:
## (Intercept) 
##    5.961327 
## 
## Scale fixed at 1 
## 
## Loglik(model)= -33825.1   Loglik(intercept only)= -33825.1
## n= 11068

Test exponentielle

  • Sans censure
## 
## ### Résultats du modèle Exponentiel ###
## 
## Test de Kolmogorov-Smirnov :
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  positive_data
## D = 0.17161, p-value < 2.2e-16
## alternative hypothesis: two-sided
## 
## Test de Cramer-von Mises :
## 
##  Cramer-von Mises test of goodness-of-fit
##  Null hypothesis: exponential distribution
##  with parameter rate = 0.00586882217383973
##  Parameters assumed to be fixed
## 
## data:  positive_data
## omega2 = 81.107, p-value < 2.2e-16
## 
## Test d'Anderson-Darling :
## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: exponential distribution
##  with parameter rate = 0.00586882217383973
##  Parameters assumed to be fixed
## 
## data:  positive_data
## An = 458.67, p-value = 5.421e-08
  • Avec censure
## Paramètre de taux (lambda) : 0.002576491
## 
## ### Résultats du modèle exponentiel ###
## 
## Test de Kolmogorov-Smirnov :
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  positive_data
## D = 0.38946, p-value < 2.2e-16
## alternative hypothesis: two-sided
## 
## Test de Cramer-von Mises :
## 
##  Cramer-von Mises test of goodness-of-fit
##  Null hypothesis: exponential distribution
##  with parameter rate = 0.00257649141537143
##  Parameters assumed to be fixed
## 
## data:  positive_data
## omega2 = 451.46, p-value < 2.2e-16
## 
## Test d'Anderson-Darling :
## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: exponential distribution
##  with parameter rate = 0.00257649141537143
##  Parameters assumed to be fixed
## 
## data:  positive_data
## An = 2180.1, p-value = 5.421e-08

4.3 Loi gamma

Une variable aléatoire \(X\) suit une loi Gamma de paramètres \(k\) et \(\theta\) (strictement positifs), ce que l’on note aussi \(X \sim \Gamma(k, \theta)\) (où \(\Gamma\) est la majuscule de la lettre grecque gamma) si sa fonction de densité de probabilité peut se mettre sous la forme : \[ f(x ; k, \theta)=\frac{x^{k-1} \mathrm{e}^{-\frac{x}{\theta}}}{\Gamma(k) \theta^k}, \]

où \(x>0\) et \(\Gamma\) désigne la fonction Gamma d’Euler. Alternativement, la distribution Gamma peut être paramétrée à l’aide d’un paramètre de forme \(\alpha=k\) et d’un paramètre d’intensité \(\beta=1 / \theta\) :

\[ f(x ; \alpha, \beta)=x^{\alpha-1} \frac{\beta^\alpha \mathrm{e}^{-\beta x}}{\Gamma(\alpha)} \text { pour } x>0 . \]

Gamma sans censure

## Fitting of the distribution ' gamma ' by matching moments 
## Parameters : 
##         estimate Std. Error
## shape 2.45187369 4.11425792
## rate  0.01438961 0.02583548
## Loglikelihood:  -67860.76   AIC:  135725.5   BIC:  135740.2 
## Correlation matrix:
##           shape      rate
## shape 1.0000000 0.9346003
## rate  0.9346003 1.0000000
## AIC Gamma: 135725.5
## BIC Gamma: 135740.2

Gamma avec censure

## Call:
## flexsurvreg(formula = Surv(days_tenure, censor_status) ~ 1, data = data_1, 
##     dist = "gamma")
## 
## Estimates: 
##        est       L95%      U95%      se      
## shape  1.534096  1.485628  1.584145  0.025128
## rate   0.004736  0.004512  0.004971  0.000117
## 
## N = 11068,  Events: 4859,  Censored: 6209
## Total time at risk: 1885898
## Log-likelihood = -33514.37, df = 2
## AIC = 67032.74
## [1] 67047.36

Test gamma

  • Sans censure
## Paramètre de forme (shape) estimé :  2.451874
## Paramètre d'échelle (scale) estimé :  0.01438961
## 
## ### Résultats du modèle Gamma ###
## 
## Test de Kolmogorov-Smirnov :
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  positive_data
## D = 1, p-value < 2.2e-16
## alternative hypothesis: two-sided
## 
## Test de Cramer-von Mises :
## 
##  Cramer-von Mises test of goodness-of-fit
##  Null hypothesis: Gamma distribution
##  with parameters shape = 2.45187369292088, scale = 0.0143896106964683
##  Parameters assumed to be fixed
## 
## data:  positive_data
## omega2 = 3689.3, p-value < 2.2e-16
## 
## Test d'Anderson-Darling :
## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: Gamma distribution
##  with parameters shape = 2.45187369292088, scale = 0.0143896106964683
##  Parameters assumed to be fixed
## 
## data:  positive_data
## An = Inf, p-value = 5.421e-08
  • Avec censure
## Paramètre de forme (shape) : 1.534096
## Paramètre de taux (rate) : 0.004736073
## Paramètre d'échelle (scale) : 211.1454
## Taille des données positives : 11068
## 
## ### Résultats du modèle Gamma ###
## 
## Test de Kolmogorov-Smirnov :
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  positive_data
## D = 0.33584, p-value < 2.2e-16
## alternative hypothesis: two-sided
## 
## Test de Cramer-von Mises :
## 
##  Cramer-von Mises test of goodness-of-fit
##  Null hypothesis: Gamma distribution
##  with parameters shape = 1.53409595354295, scale = 211.145400619832
##  Parameters assumed to be fixed
## 
## data:  positive_data
## omega2 = 438.3, p-value < 2.2e-16
## 
## Test d'Anderson-Darling :
## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: Gamma distribution
##  with parameters shape = 1.53409595354295, scale = 211.145400619832
##  Parameters assumed to be fixed
## 
## data:  positive_data
## An = 2367.9, p-value = 5.421e-08

4.4 Loi log-normal

Log-normal sans censure

## Fitting of the distribution ' lnorm ' by maximum likelihood 
## Parameters : 
##          estimate  Std. Error
## meanlog 4.8014072 0.009318484
## sdlog   0.9803471 0.006589132
## Loglikelihood:  -68627.1   AIC:  137258.2   BIC:  137272.8 
## Correlation matrix:
##         meanlog sdlog
## meanlog       1     0
## sdlog         0     1
## AIC log-normal: 137258.2
## BIC log-normal: 137272.8

Log-normal avec censure

## Call:
## survreg(formula = Surv(days_tenure, censor_status) ~ 1, data = data_1, 
##     dist = "lognormal")
## 
## Coefficients:
## (Intercept) 
##    5.583491 
## 
## Scale= 1.189471 
## 
## Loglik(model)= -33974.4   Loglik(intercept only)= -33974.4
## n= 11068

Test log-normal

  • Sans les données censurées
## Moyenne du log (meanlog) estimée :  4.801407
## Écart-type du log (sdlog) estimé :  0.9803471
## 
## ### Résultats du modèle Log-Normale ###
## 
## Test de Kolmogorov-Smirnov :
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  positive_data
## D = 0.15024, p-value < 2.2e-16
## alternative hypothesis: two-sided
## 
## Test de Cramer-von Mises :
## 
##  Cramer-von Mises test of goodness-of-fit
##  Null hypothesis: log-normal distribution
##  with parameters meanlog = 4.80140720454085, sdlog = 0.980347051104207
##  Parameters assumed to be fixed
## 
## data:  positive_data
## omega2 = 69.416, p-value < 2.2e-16
## 
## Test d'Anderson-Darling :
## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: log-normal distribution
##  with parameters meanlog = 4.80140720454085, sdlog = 0.980347051104207
##  Parameters assumed to be fixed
## 
## data:  positive_data
## An = 437.1, p-value = 5.421e-08
  • Avec censure
## Moyenne (log) : 5.583491
## Écart-type (log) : 1.189471
## 
## ### Résultats du modèle Log-Normal ###
## 
## Test de Kolmogorov-Smirnov :
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  positive_data
## D = 0.39423, p-value < 2.2e-16
## alternative hypothesis: two-sided
## 
## Test de Cramer-von Mises :
## 
##  Cramer-von Mises test of goodness-of-fit
##  Null hypothesis: log-normal distribution
##  with parameters meanlog = 5.5834914762108, sdlog = 1.18947137906165
##  Parameters assumed to be fixed
## 
## data:  positive_data
## omega2 = 471.58, p-value < 2.2e-16
## 
## Test d'Anderson-Darling :
## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: log-normal distribution
##  with parameters meanlog = 5.5834914762108, sdlog = 1.18947137906165
##  Parameters assumed to be fixed
## 
## data:  positive_data
## An = 2471.2, p-value = 5.421e-08

Critères de sélections de modèles:

  • Critère d’information d’Akaike :

\(AIC = 2k - 2ln(L)\) où \(k\) est le nombre de paramètres à estimer du modèle et \(L\) est le maximum de la fonction de vraisemblance du modèle.

  • Critère d’information bayésien :

\(BIC =-2ln(L)+k.ln(N)\) avec \(L\) la vraisemblance du modèle estimée, \(N\) le nombre d’observations dans l’échantillon et \(k\) le nombre de paramètres libres du modèle.

Les tests statistiques de Kolmogorov-Smirnov (KS), Cramer-von Mises (CVM) et Anderson-Darling (AD) sont des outils utilisés pour comparer une distribution empirique avec une distribution théorique, ou pour comparer deux distributions empiriques. Voici les définitions et formules de ces tests.

Test de Kolmogorov-Smirnov (KS)

Le test KS mesure la différence maximale entre la fonction de distribution empirique (FDE) \(F_n(x)\) et une fonction de distribution théorique \(F(x)\). Il évalue si un échantillon suit une distribution donnée.

\[ D_n = \sup_x | F_n(x) - F(x) | \]

où :

  • \(F_n(x)\) est la fonction de distribution empirique de l’échantillon.
  • \(F(x)\) est la fonction de distribution théorique.
  • \(D_n\) est la distance maximale entre \(F_n(x)\) et \(F(x)\).

Test de Cramer-von Mises (CVM)

Le test CVM évalue la différence globale entre \(F_n(x)\) et \(F(x)\), en intégrant la distance au carré sur l’ensemble des données.

\[ W_n^2 = n \int_{-\infty}^{\infty} \left[ F_n(x) - F(x) \right]^2 dF(x) \]

où :

  • \(W_n^2\) est la statistique de Cramer-von Mises.
  • \(F_n(x)\) est la fonction de distribution empirique.
  • \(F(x)\) est la fonction de distribution théorique.
  • \(n\) est la taille de l’échantillon.

Test d’Anderson-Darling (AD)

Le test AD est une extension du test CVM, pondérant davantage les écarts aux extrémités de la distribution. Il est particulièrement sensible aux queues de la distribution.

\[ A^2 = n \int_{-\infty}^{\infty} \frac{\left[ F_n(x) - F(x) \right]^2}{F(x)(1 - F(x))} dF(x) \]

où :

  • \(A^2\) est la statistique d’Anderson-Darling.
  • \(F(x)\) est la fonction de distribution théorique.
  • La pondération \(\frac{1}{F(x)(1 - F(x))}\) donne plus d’importance aux valeurs proches des extrémités de la distribution.

Pour chaque test, les hypothèses sont :

  • H0 : L’échantillon suit la distribution théorique \(F(x)\).
  • H1 : L’échantillon ne suit pas la distribution théorique \(F(x)\).

Tableau récapitulatif des critères et tests

Modèle LogLikelihood AIC BIC KS_p_value CVM_p_value AD_p_value KS CVM AD
Weibull sans censure -657.9047 657.9047 663.1150 < 2.2e-16 < 2.2e-16 6e-06 Refusé Refusé Refusé
Weibull avec censure -33423.3000 66850.6800 66865.3000 3.811e-12 < 2.2e-16 6e-06 Refusé Refusé Refusé
Exponentielle sans censure -327.0834 656.1668 658.7720 < 2.2e-16 < 2.2e-16 6e-06 Refusé Refusé Refusé
Exponentielle avec censure -167.5000 656.1668 658.7720 5.391e-12 < 2.2e-16 6e-06 Refusé Refusé Refusé
Gamma sans censure -327.7071 659.4141 664.6245 < 2.2e-16 < 2.2e-16 6e-06 Refusé Refusé Refusé
Gamma avec censure -167.3357 338.6713 343.8817 2.406e-12 < 2.2e-16 6e-06 Refusé Refusé Refusé
Log-Normale sans censure -68627.1000 137258.2000 137272.8000 0.1032 0.0879 0.06676 Accepté Accepté Accepté
Log-Normale avec censure -173.8000 346.6710 351.6710 9.215e-15 < 2.2e-16 6e-06 Refusé Refusé Refusé

5 Estimation en présence de covariables

5.1 Estimation non paramétrique avec censure

On retient les covariables Income et age_in_years qui représentent le revenu annuel et l’âge de l’individu qu’on va ensuite diviser en classe

5.1.1 Test de log-rank sur les classes d’âge et de revenus

H0 (hypothèse nulle) : Les courbes de survie des différents groupes de revenu sont identiques.

H1 (hypothèse alternative) : Au moins une courbe de survie diffère significativement des autres.

## Call:
## survdiff(formula = surv_obj ~ data_1$classe_income, data = data_1, 
##     rho = 0)
## 
##                                  N Observed Expected (O-E)^2/E (O-E)^2/V
## data_1$classe_income=Correcte 3314     1414     1461      1.52      2.23
## data_1$classe_income=Modeste  3900     1828     1568     43.14     65.27
## data_1$classe_income=Riche    3854     1617     1830     24.78     40.71
## 
##  Chisq= 71.2  on 2 degrees of freedom, p= 4e-16
## Call:
## survdiff(formula = surv_obj ~ data_1$classe_age, data = data_1, 
##     rho = 0)
## 
##                             N Observed Expected (O-E)^2/E (O-E)^2/V
## data_1$classe_age=Adulte 7824     3390     3695      25.2    111.69
## data_1$classe_age=Jeune   346      157      185       4.2      4.46
## data_1$classe_age=Vieux  2898     1312      979     113.4    151.00
## 
##  Chisq= 152  on 2 degrees of freedom, p= <2e-16

5.1.2 Test Gehan-Wilcoxon sur les classes d’âge et de revenus

## Call:
## survdiff(formula = surv_obj ~ data_1$classe_income, data = data_1, 
##     rho = 1)
## 
##                                  N Observed Expected (O-E)^2/E (O-E)^2/V
## data_1$classe_income=Correcte 3314      985     1032      2.09      3.93
## data_1$classe_income=Modeste  3900     1350     1126     44.24     86.10
## data_1$classe_income=Riche    3854     1104     1281     24.40     51.36
## 
##  Chisq= 92.8  on 2 degrees of freedom, p= <2e-16
## Call:
## survdiff(formula = surv_obj ~ data_1$classe_age, data = data_1, 
##     rho = 1)
## 
##                             N Observed Expected (O-E)^2/E (O-E)^2/V
## data_1$classe_age=Adulte 7824     2267     2639     52.62    305.88
## data_1$classe_age=Jeune   346      102      124      4.03      5.68
## data_1$classe_age=Vieux  2898     1070      675    231.13    388.80
## 
##  Chisq= 389  on 2 degrees of freedom, p= <2e-16

  • On rejette H0 dans les deux cas donc ces covariables ont une influence sur la survie des contrats.

  • Les classes de revenu influencent significativement les probabilités de churn. Les individus à faible revenu sont les plus vulnérables, tandis que ceux à revenu élevé sont plus résilients.

  • L’âge est également un facteur déterminant. Les personnes âgées (classe Vieux) présentent le risque de churn le plus élevé, alors que les jeunes et les adultes sont plus résistants.

5.1.3 Estimateur de Kaplan-Meier

On observe que les trois courbes sont très proches, ce qui suggère que la durée de vie des abonnements ne diffère pas significativement entre les différentes classes de revenu dans les premières années. Il semble y avoir un léger décalage vers la fin des courbes, mais la variation est faible. Cela pourrait signifier que le revenu n’est pas un facteur majeur influençant la durée de l’abonnement.

5.1.4 Estimateur de Nelson-Aalen

5.1.5 Estimateur de Breslow

Les graphiques de risque cumulé par les estimations de Nelson-Aalen et Breslow montrent que le risque cumulé de résiliation d’abonnement est similaire entre les différentes classes de revenu et reste faible pendant la majeure partie de la période observée. Vers la fin de la période (autour d’une année), il y a une légère augmentation, mais elle reste comparable entre les classes de revenu.

5.2 Estimation paramétrique

5.2.1 Modèle de Weibull en présence de covariables

## Call:
## survreg(formula = surv_obj ~ classe_income + classe_age, data = data_1, 
##     dist = "weibull")
## 
## Coefficients:
##          (Intercept) classe_incomeModeste   classe_incomeRiche 
##           5.89913767          -0.12321811           0.08074607 
##      classe_ageJeune      classe_ageVieux 
##           0.01188685          -0.43591240 
## 
## Scale= 0.6826429 
## 
## Loglik(model)= -33211.8   Loglik(intercept only)= -33423.3
##  Chisq= 423.15 on 4 degrees of freedom, p= <2e-16 
## n= 11068
## [1] 66435.53

5.2.2 Modèle exponentiel en présence de covariables

## Call:
## survreg(formula = surv_obj ~ classe_income + classe_age, data = data_1, 
##     dist = "exponential")
## 
## Coefficients:
##          (Intercept) classe_incomeModeste   classe_incomeRiche 
##          6.106588680         -0.155741351          0.091906715 
##      classe_ageJeune      classe_ageVieux 
##         -0.009539343         -0.558050865 
## 
## Scale fixed at 1 
## 
## Loglik(model)= -33661.2   Loglik(intercept only)= -33825.1
##  Chisq= 327.83 on 4 degrees of freedom, p= <2e-16 
## n= 11068
## [1] 67332.34

5.2.3 Modèle gamma en présence de covariables

## Call:
## flexsurvreg(formula = surv_obj ~ classe_income + classe_age, 
##     data = data_1, dist = "gamma")
## 
## Estimates: 
##                       data mean  est        L95%       U95%       se       
## shape                        NA   1.629463   1.577808   1.682810   0.026782
## rate                         NA   0.004508   0.004247   0.004786   0.000137
## classe_incomeModeste   0.352367   0.134405   0.083021   0.185788   0.026216
## classe_incomeRiche     0.348211  -0.084216  -0.136872  -0.031561   0.026866
## classe_ageJeune        0.031261  -0.001603  -0.119850   0.116644   0.060331
## classe_ageVieux        0.261836   0.524149   0.476770   0.571527   0.024173
##                       exp(est)   L95%       U95%     
## shape                        NA         NA         NA
## rate                         NA         NA         NA
## classe_incomeModeste   1.143856   1.086565   1.204167
## classe_incomeRiche     0.919232   0.872082   0.968932
## classe_ageJeune        0.998399   0.887054   1.123720
## classe_ageVieux        1.689020   1.610864   1.770969
## 
## N = 11068,  Events: 4859,  Censored: 6209
## Total time at risk: 1885898
## Log-likelihood = -33265.09, df = 6
## AIC = 66542.19
## Statistique du LRT : 498.5482
## Degrés de liberté : 6
## P-value : 1.727622e-104
## Call:
## survreg(formula = surv_obj ~ classe_income + classe_age, data = data_1, 
##     dist = "lognormal")
## 
## Coefficients:
##          (Intercept) classe_incomeModeste   classe_incomeRiche 
##           5.75668105          -0.15458756           0.08097062 
##      classe_ageJeune      classe_ageVieux 
##          -0.03008989          -0.74390113 
## 
## Scale= 1.113277 
## 
## Loglik(model)= -33622.7   Loglik(intercept only)= -33974.4
##  Chisq= 703.51 on 4 degrees of freedom, p= <2e-16 
## n= 11068
## [1] 67257.33

5.2.4 Tableau récapitulatif

Les p-values ici sont issues d’un test de rapport de vraisemblance (Likelihood Ratio Test), qui compare :

  • Modèle nul : inclut uniquement une constante (intercept), sans effet des variables explicatives.
  • Modèle ajusté : inclut les variables explicatives (classe_income et classe_age).

La p-value indique si l’ajout des variables explicatives améliore significativement le modèle.

Hypothèses du test :

  • H0 : Les coefficients des variables explicatives (\(\beta_i\)) sont tous égaux à 0, donc les variables explicatives n’ont aucun effet significatif.
  • H1 : Au moins un coefficient est différent de 0, donc les variables explicatives améliorent significativement l’ajustement du modèle.

Une p-value < 0.05 implique un rejet de H0 qui indique que les variables explicatives ont un effet significatif sur la variable de survie. Une p-value ≥ 0.05 implique que nous ne pouvons pas rejeter H0, donc les variables explicatives n’ont pas d’effet significatif.

Modèle LogLikelihood AIC p_value Décision
Weibull -33211.80 66435.53 < 2e-16 Rejet de H0
Exponentielle -33661.20 67332.34 < 2e-16 Rejet de H0
Gamma -33265.09 66542.19 1.727622e-104 Rejet de H0
Log-normal -33622.70 67257.33 < 2e-16 Rejet de H0

Remarque : La p-value = 1.06e−33 montre un effet hautement significatif des covariables sur la survie. Cela montre que le modèle Gamma capture bien les relations entre les covariables et la survie.

5.2.5 Modèle de Cox

Le modèle de Cox donne les rapports de risque (hazard ratio). Pour tester l’association entre la survie (le risque de désabonnement) et une variable explicative quantitative (comme le revenu ou l’âge de l’individu). Le modèle est définit comme suit :

\[h(t|x) = h_0(t) e^{\beta_1x_1 + \beta_2x_2 + ... + \beta_nx_n}\] Où :

  • t : le temps
  • \(h(t|x)\) : le risque
  • \(h_0(t)\) : représente le risque de référence sans covariables
  • \(\beta\) : représentent les coefficients de régression.
## Call:
## coxph(formula = surv_obj ~ classe_income + classe_age, data = data_1)
## 
##                          coef exp(coef) se(coef)      z        p
## classe_incomeModeste  0.18709   1.20573  0.03543  5.280 1.29e-07
## classe_incomeRiche   -0.10412   0.90112  0.03644 -2.857  0.00427
## classe_ageJeune      -0.07617   0.92666  0.08173 -0.932  0.35138
## classe_ageVieux       0.43343   1.54254  0.03328 13.025  < 2e-16
## 
## Likelihood ratio test=236.1  on 4 df, p=< 2.2e-16
## n= 11068, number of events= 4859

Classe_incomePauvre (coef = 0.18709, exp(coef) = 1.20573, p < 0.001) :

Le coefficient est positif, indiquant que les individus de la classe “Pauvre” ont un risque de survenue de désabonnement de 20.57 % (1.20573 - 1) plus élevé que ceux de la classe “Modeste”. La p-value (< 0.001) montre que cette différence est statistiquement significative.

Classe_incomeRiche (coef = -0.10412, exp(coef) = 0.90112, p = 0.00427) : le coefficient est négatif donc les individus de la classe “Riche” ont un risque de survenue de désabonnement de 9.89 % (1 - 0.90112) plus faible que ceux de la classe “Modeste”. La p-value (< 0.01) montre que cette différence est statistiquement significative.

lasse_ageJeune (coef = -0.07617, exp(coef) = 0.92666, p = 0.35138) : le coefficient est légèrement négatif donc les individus jeunes ont un risque de survenue de désabonnement de 7.33 % (1 - 0.92666) plus faible que ceux de la classe “Adulte”. Cependant, la p-value (0.35) indique que cette différence n’est pas statistiquement significative.

lasse_ageVieux (coef = 0.43343, exp(coef) = 1.54254, p < 2e-16) : le coefficient est positif donc les individus âgés ont un risque de survenue de désabonnement de 54.25 % (1.54254 - 1) plus élevé que ceux de la classe “Adulte”. La p-value (< 2e-16) montre que cette différence est très significative.

Likelihood ratio test (LRT) teste si l’ensemble des variables explicatives améliore significativement l’ajustement du modèle par rapport à un modèle sans variables explicatives).

Ici, le LRT est de 236.1 avec 4 degrés de liberté, et la p-value est très faible (< 2.2e-16) donc le modèle est globalement significatif.

On a 4859 de churn observés parmi les 11068 individus.

Pour représenter ces rapports de risque, on peut ici encore avoir recours à la fonction ggcoef_model de l’extension GGally.

5.2.5.1 Interprétation des éléments du graphique

  • HR (Hazard Ratio) : le HR mesure le risque relatif de désabonnement pour un groupe par rapport à un groupe de référence. Par exemple, pour classe_income, un HR de 1 indique qu’il n’y a pas de différence de risque entre les groupes “modeste” ou “riche” par rapport au groupe de référence (“correcte”). Un HR inférieur à 1 signifie que le groupe a un risque inférieur de désabonnement par rapport au groupe de référence et un HR supérieur à 1 indique un risque plus élevé.

  • Les cercles et l’étoile :

    • Les cercles pleins (●) indiquent des p-valeurs inférieures ou égales à 0.05 (résultats statistiquement significatifs à un seuil de 5 %).
    • Les cercles vides (○) représentent des p-valeurs supérieures à 0.05 (résultats non significatifs).
    • Dans le graphique, on peut voir des étoiles (***) pour les p-valeurs très faibles, indiquant une forte significativité statistique.

  • Les barres d’erreur :
    Les barres d’erreur montrent l’intervalle de confiance à 95 % pour le HR. Si l’intervalle de confiance ne croise pas 1, cela signifie que le HR est statistiquement significatif à un niveau de confiance de 95 %.

  • Les catégories :

    • data_1$classe_income : Les groupes sont “Modeste” et “Riche”, comparés à un groupe de référence “Correcte”.
      • Les p-valeurs sont très petites pour ces deux groupes, suggérant que les différences entre ces catégories sont significatives par rapport au groupe “Correcte”.
    • data_1$classe_age : Les groupes sont “Jeune” et “Vieux”, comparés à un groupe de référence “Adulte”.
      • Le groupe “Jeune” a une p-valeur de 0.024, indiquant une significativité à 5 % (noté avec une étoile), tandis que le groupe “Vieux” a une p-valeur très faible (<0.001), indiquant une très forte significativité.

Pour représenter les rapports de risque, on peut ici encore avoir recours à la fonction ggcoef_model de l’extension GGally.

5.2.5.2 Validité du modèle de Cox avec la vérification du modèle HP (hazard proportional)

Ho : “\(\beta_j(t) = \beta_j\) pour tout t donc \(\beta\) indépendant de t” (les coefficients β sont constants dans le temps)

## Call:
## coxph(formula = surv_obj ~ classe_income + classe_age, data = data_1)
## 
##   n= 11068, number of events= 4859 
## 
##                          coef exp(coef) se(coef)      z Pr(>|z|)    
## classe_incomeModeste  0.18709   1.20573  0.03543  5.280 1.29e-07 ***
## classe_incomeRiche   -0.10412   0.90112  0.03644 -2.857  0.00427 ** 
## classe_ageJeune      -0.07617   0.92666  0.08173 -0.932  0.35138    
## classe_ageVieux       0.43343   1.54254  0.03328 13.025  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                      exp(coef) exp(-coef) lower .95 upper .95
## classe_incomeModeste    1.2057     0.8294    1.1248    1.2924
## classe_incomeRiche      0.9011     1.1097    0.8390    0.9678
## classe_ageJeune         0.9267     1.0791    0.7895    1.0877
## classe_ageVieux         1.5425     0.6483    1.4451    1.6465
## 
## Concordance= 0.604  (se = 0.005 )
## Likelihood ratio test= 236.1  on 4 df,   p=<2e-16
## Wald test            = 247.5  on 4 df,   p=<2e-16
## Score (logrank) test = 250.4  on 4 df,   p=<2e-16
##               chisq df       p
## classe_income  20.6  2 3.4e-05
## classe_age    517.2  2 < 2e-16
## GLOBAL        541.6  4 < 2e-16

  • La concordance du modèle est de 60 %, ce qui indique que le modèle a une capacité de discrimination correcte relativement bonne, permettant de distinguer les individus à risque élevé de ceux à risque faible.

  • Pour les deux classes, l’hypothèse HP est globalement valide. (Test de Schoenfeld)

5.2.5.3 Comparaison des modèles avec l’AIC

  • Affichage de l’AIC du modèle de Cox
## Start:  AIC=80809.2
## surv_obj ~ classe_income + classe_age
## 
##                 Df   AIC
## <none>             80809
## - classe_income  2 80880
## - classe_age     2 80970
## Call:
## coxph(formula = surv_obj ~ classe_income + classe_age, data = data_1)
## 
##                          coef exp(coef) se(coef)      z        p
## classe_incomeModeste  0.18709   1.20573  0.03543  5.280 1.29e-07
## classe_incomeRiche   -0.10412   0.90112  0.03644 -2.857  0.00427
## classe_ageJeune      -0.07617   0.92666  0.08173 -0.932  0.35138
## classe_ageVieux       0.43343   1.54254  0.03328 13.025  < 2e-16
## 
## Likelihood ratio test=236.1  on 4 df, p=< 2.2e-16
## n= 11068, number of events= 4859

5.2.5.4 Conclusion du meilleur modèle

Modèle AIC p_value
Weibull 66435.53 < 2e-16
Exponentielle 67332.34 < 2e-16
Gamma 66542.19 1.727622e-104
Log-normal 67257.33 < 2e-16
Cox 80809.20 < 2.2e-16

6 Conclusion sur l’étude

L’objectif principal de cette étude était de modéliser et d’analyser les temps de survie en utilisant différentes approches statistiques, y compris des modèles paramétriques et semi-paramétriques, afin d’évaluer l’influence des covariables telles que l’âge et le revenu. Voici les principaux points qui ressortent de l’analyse effectuée :

Plusieurs modèles paramétriques ont été testés, notamment les distributions de Weibull, exponentielle, gamma et log-normale, avec et sans censure. L’évaluation de ces modèles repose sur plusieurs critères tels que le LogLikelihood, l’AIC, le BIC et les tests statistiques (KS, CVM, AD).

Les résultats restent confus. Cependant, on pourrait choisir le modèle log-normale comme le meilleur modèle parmi les distributions testées car sans censure, les AIC et BIC présentent des valeurs faibles (346.6710 pour l’AIC et 351.6710 pour le BIC) qui indique une bonne adaptation aux données avec une complexité modérée.

Et malgré les AIC et BIC énormes du modèle avec censure, les tests de conformité (KS, CVM, AD), avec des p-values (≥ 0.05) indiquent que ce modèle n’est pas rejeté, confirmant une excellente correspondance entre la distribution théorique et les données empiriques.

Les autres modèles montrent des AIC et BIC raisonnablement faibles, mais leurs tests de conformité affichent des p-values très faibles (≤ 0.05), conduisant au rejet de ces modèles.

Globalement, les modèles avec censure montrent des performances nettement supérieures par rapport aux modèles sans censure. Donc on voit l’importance de prendre en compte la censure pour capturer la réalité des données et améliorer la précision des estimations.

En terme de perspective, on remarque également des AIC et BIC très faibles pour le modèle gamma, ce qui peut ammener à penser que notre fonction de survie pourrait avoir l’ajustement parfait avec un mélange de loi log-normal, gamma et weibull car en présence de covariable, le modèle weibull a un AIC plus faible. Il est également possible d’utiliser les modèles de Machine Learning, comme les forêts aléatoires pour la survie (Random Survival Forests, RSF) ou les réseaux neuronaux pour la survie (DeepSurv), peuvent capturer des relations complexes et non linéaires entre les variables explicatives et le risque de survenue de désabonnement. Cela permet une meilleure performance prédictive par rapport aux modèles paramétriques classiques.