Sujet 6: Autour du Paradoxe de Simpson
+Table des matières
+ +1 Introduction
++En 1972-1974, à Whickham, une ville du nord-est de l'Angleterre, située à environ 6,5 kilomètres au sud-ouest de Newcastle upon Tyne, un sondage d'un sixième des électeurs a été effectué afin d'éclairer des travaux sur les maladies thyroïdiennes et cardiaques (Tunbridge et al. 1977). Une suite de cette étude a été menée vingt ans plus tard (Vanderpump et al. 1995). Certains des résultats avaient trait au tabagisme et cherchaient à savoir si les individus étaient toujours en vie lors de la seconde étude. Par simplicité, nous nous restreindrons aux femmes et parmi celles-ci aux 1314 qui ont été catégorisées comme "fumant actuellement" ou "n'ayant jamais fumé". Il y avait relativement peu de femmes dans le sondage initial ayant fumé et ayant arrêté depuis (162) et très peu pour lesquelles l'information n'était pas disponible (18). La survie à 20 ans a été déterminée pour l'ensemble des femmes du premier sondage. +
+ ++Vous trouverez sur chaque ligne si la personne fume ou non, si elle est vivante ou décédée au moment de la seconde étude, et son âge lors du premier sondage. +(source) +
+2 Mission
+2.1 Relation entre tabagisme et décès
++Représentez dans un tableau le nombre total de femmes vivantes +et décédées sur la période en fonction de leur habitude de +tabagisme. Calculez dans chaque groupe (fumeuses / non +fumeuses) le taux de mortalité (le rapport entre le nombre de +femmes décédées dans un groupe et le nombre total de femmes +dans ce groupe). Vous pourrez proposer une représentation +graphique de ces données et calculer des intervalles de +confiance si vous le souhaitez. En quoi ce résultat est-il +surprenant ? +
+ ++Tout d'abord Chargons le jeu de données dans une dataframe pandas pour +faciliter la manipulation de données +
+from pathlib import Path +import matplotlib.pyplot as plt +import pandas as pd + +# set floating precision to 1 (ie: 97.244->97.2) +pd.set_option("display.precision", 1) + +data_url = "http://gitlab.inria.fr/learninglab/mooc-rr/mooc-rr-ressources/blob/master/module3/Practical_session/Subject6_smoking.csv" +data_file = Path("Subject6_smoking.csv") + +if not data_file.exists(): + raise IOError(f"Jeu de données introuvable, merci de le télécharger à nouveau via ce lien {data_url}") + +data = pd.read_csv(data_file, dtype={'Smoker': 'category', 'Status': 'category'}) +data ++
+ Smoker Status Age +0 Yes Alive 21.0 +1 Yes Alive 19.3 +2 No Dead 57.5 +3 No Alive 47.1 +4 Yes Alive 81.4 +... ... ... ... +1309 Yes Alive 35.9 +1310 No Alive 22.3 +1311 Yes Dead 62.1 +1312 No Dead 88.6 +1313 No Alive 39.1 + +[1314 rows x 3 columns] ++ +
data.describe(include='all')
+
++ Smoker Status Age + count 1314 1314 1314.0 + unique 2 2 NaN + top No Alive NaN + freq 732 945 NaN + mean NaN NaN 47.4 + std NaN NaN 19.2 + min NaN NaN 18.0 + 25% NaN NaN 31.3 + 50% NaN NaN 44.8 + 75% NaN NaN 60.6 + max NaN NaN 89.9 ++ +
+On a bien deux valeurs possibles pour les variables catégorielles: +
+-
+
- Fumeuse(Smoker): Yes/No
+
-
+
- Indique si la personne est fumeuse au moment de la seconde étude +(20 ans après la première) +
+ - Etat(Status) : Alive/Dead
+
-
+
- Indique si la personne est vivante au moment de la seconde étude +
+
+Au moment de la première étude, les 1314 femmes sollicitées étaient majeures. La moitié avait moins de 44,8 ans et +la doyenne était agée de 89,9 ans. +
+ ++Etudions le lien entre tabagisme et décès via ce tableau de +contingence sur les variables Smoker et Status: +
+pd.crosstab(data.Smoker, data.Status, margins=True)
+
++Status Alive Dead All +Smoker +No 502 230 732 +Yes 443 139 582 +All 945 369 1314 + ++ +
+Depuis la première étude, 369 personnes sont décédées(\(39\%=\frac{369}{1314}\)). On observe +aussi que le nombre de non fumeuses est sensiblement plus élévé(\(+26\%=\frac{732}{582}\)) +que le nombre de fumeuses. Ainsi, on peut s'attendre à avoir une +distortion sur le nombre de morts. C'est effectivement le cas, mais +bien plus qu'on ne le supposerait car les non-fumeuses sont décédées +en bien plus grand nombre (\(+65%\%=\frac{230}{139}\)) que les fumeuses! +
+ ++Regardons maintenant les mêmes données sous forme de ratios par +rapport à toute la population. Par exemple, 10,6% des femmes de +l'étude étaient fumeuses avant de décéder. +
+pd.crosstab(data.Smoker, data.Status, margins=True, normalize='all')*100 ++
+Status Alive Dead All +Smoker +No 38.2 17.5 55.7 +Yes 33.7 10.6 44.3 +All 71.9 28.1 100.0 + ++ +
+Les ratios confirment ce que nous avons vu précédemment: les femmes +sont inégalement réparties selon la modalité fumeuse(44,3%)/non fumeuse(55,7%) +
+ ++Déterminons maintenant les taux de mortalités pour chaque groupe +(fumeuse/ non fumeuse): +
+pd.crosstab(data.Smoker, data.Status, margins=True, normalize='index')*100 ++
+Status Alive Dead +Smoker +No 68.6 31.4 +Yes 76.1 23.9 +All 71.9 28.1 + ++ +
+Contre-intuitivement, nous observons que proportionnellement le taux +de mortalité est bien plus élevé chez les non-fumeuses que chez les +fumeuses(31,4% contre 23,9%) +
+ ++Observons ce résultat surprenant sous forme graphique: +
+ +pd.crosstab(data.Smoker, data.Status).plot.barh(stacked=True) + +figure_path = "data_as_stacked_bar_plot.png" +plt.savefig(figure_path) +figure_path ++
+
+Qu'est-ce qui pourrait justifier que la surmortalité chez les +non-fumeuses ? +Existerait-il une autre variable qui +influencerait ces résultats ? Une autre comorbidité par exemple ? +Ou alors fumeraient-elle des MELL-O-WELL comme le recommande ce +médecin sur une publicité de 1930 ?! +
+ +
+
+(source)
+
2.2 Partitionnement selon la classe d'âge
++Reprenez la question 1 (effectifs et taux de mortalité) +en rajoutant une nouvelle catégorie liée à la classe d'âge. On +considérera par exemple les classes suivantes : 18-34 ans, +34-54ans, 55-64 ans, plus de 65 ans.En quoi ce résultat est-il +surprenant ? Arrivez-vous à expliquer ceparadoxe ? De même, +vous pourrez proposer une représentationgraphique de ces +données pour étayer vos explications. +
+ + ++Comme suggéré, groupons par tranches d'age les femmes de l'études: +
+# On choisit 17 ans en borne inférieur pour être sûr d'inclure le minimum (18) +data['AgeCut'] = pd.cut(data.Age, bins=[17, 34, 54, 64, data.Age.max()], labels=["18-34 ans", "34-54 ans", "55-64 ans", "+ de 65 ans"]) +data.AgeCut ++
+0 18-34 ans +1 18-34 ans +2 55-64 ans +3 34-54 ans +4 + de 65 ans + ... +1309 34-54 ans +1310 18-34 ans +1311 55-64 ans +1312 + de 65 ans +1313 34-54 ans +Name: AgeCut, Length: 1314, dtype: category +Categories (4, object): [18-34 ans < 34-54 ans < 55-64 ans < + de 65 ans] ++ +
+Remarquons que le choix des tranches d'age est ici arbitraire. On +aurait pu choisir ces 4 tranches de manière homogènes en regroupant par +quartile: +
+print("Avec Quartiles:", *pd.qcut(data.Age, 4, retbins=True)[1].astype(int), sep='|') +print("Sans Quartiles:", *[18, 34, 54, 64, 89], sep='|') ++
+Avec Quartiles:|18|31|44|60|89 +Sans Quartiles:|18|34|54|64|89 + ++ +
+Dénombrons les femmes ayant survécu/décédé selon leur tranche d'age +et si elles fumaient ou non: +
+pd.crosstab([data.Smoker, data.AgeCut], data.Status, margins=True)
+
++Status Alive Dead All +Smoker AgeCut +No 18-34 ans 213 6 219 + 34-54 ans 180 19 199 + 55-64 ans 81 40 121 + + de 65 ans 28 165 193 +Yes 18-34 ans 176 5 181 + 34-54 ans 196 41 237 + 55-64 ans 64 51 115 + + de 65 ans 7 42 49 +All 945 369 1314 ++ +
+Si on regarde maintenant les ratios, on se rend compte qu'il y avait +sensiblement moins de fumeuses de plus de 65 ans que de fumeuses (3,7% +contre 14,7%) lors de la première étude. +
+ ++Aussi, la mortalité semble plus élevée chez les fumeuses des tranches de +34 à 54 ans(3,1% contre 1,4%) +et de 55 à 64 ans(3,9% contre 3,0%). +Même si on a vu que la classe fumeuse/non fumeuse est déséquilibrée, +c'est étrange que cette tendance s'inverse sévèrement en passant +dans la tranche des plus de 65 ans(3,2% contre 12,6%) +Faudrait-il fumer toute sa vie pour enfin voir des effets positifs sur sa +santé ?! +
+ +pd.crosstab([data.AgeCut, data.Smoker], data.Status, margins=True, normalize='all')*100 ++
+Status Alive Dead All +AgeCut Smoker +18-34 ans No 16.2 0.5 16.7 + Yes 13.4 0.4 13.8 +34-54 ans No 13.7 1.4 15.1 + Yes 14.9 3.1 18.0 +55-64 ans No 6.2 3.0 9.2 + Yes 4.9 3.9 8.8 ++ de 65 ans No 2.1 12.6 14.7 + Yes 0.5 3.2 3.7 +All 71.9 28.1 100.0 ++ +
+L'observation précédente n'est pas tout à fait exacte car elle compare +les tendances sur le ratio de la population globale. Si on se concentre sur les groupes +(fumeuse/non-fumeuse, tranche d'age) on se rend compte que les taux de +mortalités n'augmentent pas si on est fumeur pour les plus jeunes(~3%) mais +aussi pour les plus agés!(~85%) +Comme précédemment, au milieu de la vie, on observe une mortalité signficativement plus +élévée chez les fumeuses: 17,3% contre 9,5% chez les 34-54 ans et +44,3% contre 33,1% chez les 55-64 ans) +Au delà de 65 ans, cette surmortalité semble disparaître. Y aurait-il +une autre variable influençant les résultats ? +
+ +pd.crosstab([data.AgeCut, data.Smoker], data.Status, margins=True, normalize='index')*100 ++
+Status Alive Dead +AgeCut Smoker +18-34 ans No 97.3 2.7 + Yes 97.2 2.8 +34-54 ans No 90.5 9.5 + Yes 82.7 17.3 +55-64 ans No 66.9 33.1 + Yes 55.7 44.3 ++ de 65 ans No 14.5 85.5 + Yes 14.3 85.7 +All 71.9 28.1 ++ +
+Visuallement, c'est plus parlant: +
+ +pd.crosstab([data.AgeCut, data.Smoker], data.Status).plot.barh(stacked=True) + +figure_path = "data_as_stacked_bar_plot_agecut.png" +plt.savefig(figure_path, bbox_inches='tight') +figure_path ++
+
+Même si le taux de mortalité semble être le même chez les femmes les plus
+agées, le nombre d'individus est considérablement plus élévé chez les
+non-fumeuses. Difficile de comparer correctement si les groupes ne sont pas homogènes.
+
+On peut supposer que la mortalité plus élevée chez les fumeuses au +cours de leur vie se poursuit au delà de 65 ans mais que d'autres +causes de mortalité viennent s'ajouter. Ce qui paraît intuitif en fin +de vie. +
+ ++Au Royaume-Uni, l'espérance de vie à la naissance en 1974 était de 72,5 ans +(source). Ainsi, lors de la première étude(de 1972 à 1974), les femmes de +plus de 65 ans pouvaient espérer vivre en moyenne 7,5ans. Ce chiffre +serait inférieur si on avait pris leur espérance de vie aux naissance +des femmes (avant 1910). +On peut naturellement attendre que 20 ans après la première étude, un +grand nombre de femmes aient décédées pour de multiples raisons, le +tabagisme n'étant que l'une des possibilités. +
+2.3 Modélisation de la probabilité de décès en fonction de l'âge
++Afin d'éviter un biais induit par des regroupements en tranches +d'âges arbitraires et non régulières, il est envisageable +d'essayer de réaliser une régression logistique. Si on +introduit une variable Death valant 1 ou 0 pour indiquer si +l'individu est décédé durant la période de 20 ans, on peut +étudier le modèle Death ~ Age pour étudier la probabilité de +décès en fonction de l'âge selon que l'on considère le groupe +des fumeuses ou des non fumeuses. Ces régressions vous +permettent-elles de conclure sur la nocivité du tabagisme? Vous +pourrez proposer une représentation graphique de ces +régressions (en n'omettant pas les régions de confiance). +
+ ++Nous avons vu précédemment que le regroupement en tranches d'age +arbitraire pouvait surréprésenter certaines modalités(fumeuse/non +fumeuse) et amener à des conclusions innattendues(les fumeuses agées +décèdent moins que les non fumeuses agées). +
+ ++En effet, nous pouvons voir sur le diagramme en violon ci-dessous que la +mortalité est plus faible chez les non fumeuses agées(+65 ans), mais +ne serait-pas parce qu'une majorité des fumeuses décèdent autour de 60 +ans ? +
+ +import seaborn as sns +plt.figure(figsize=(9, 7)) +sns.violinplot(x='Smoker', y='Age', hue='Status', data=data, split=True, inner='quartile') +sns.despine() +figure_path = "smoker_survival_violin.png" +plt.savefig(figure_path, bbox_inches='tight') +figure_path ++
+
+Comme conseillé, nous créeons une variable Death à partir de la +variable Status. Il s'agit simplement de faire correspondre 'Alive' à +0 et 'Dead' à 1: +
+# On définit une variable Death(0 si survie, 1 si décès) +data['Death'] = data.Status.cat.codes +data[['Status', 'Death']] ++
+ Status Death +0 Alive 0 +1 Alive 0 +2 Dead 1 +3 Alive 0 +4 Alive 0 +... ... ... +1309 Alive 0 +1310 Alive 0 +1311 Dead 1 +1312 Dead 1 +1313 Alive 0 + +[1314 rows x 2 columns] ++ + +
+Observons rapidement ce qu'une modélisation avec une régression +logistique sur chaque groupe fumeuse/non fumeuse pourrait donner: +
+plt.close() +plt.figure(figsize=(12, 15)) +sns.lmplot(x="Age", y="Death", data=data, logistic=True, hue='Smoker', y_jitter=.03, aspect=1.5, markers=['o', 'x']); + +sns.despine() +figure_path = "smoker_survival_logistic.png" +plt.savefig(figure_path, bbox_inches='tight') +figure_path ++
+
+Plusieurs remarques intéressantes: +
+-
+
- La modélisation en sigmoïde(logistique) chez les non-fumeuses semble +mieux coller aux observations que chez les fumeuses. Vers les 60 +ans, on observe une courbure tendant vers un point d'inflexion(ce +qui va dans le sens de la régression logistique). En revanche, dans +cette même zône, on observe une courbure linéaire chez les +fumeuses. +
- Plus l'age des fumeuses augmente, plus l'intervalle de confiance est +large. Cela révèle une fois de plus le manque de données concernant +les personnes agées chez les fumeuses +
- L'intervalle de confiance semble plus étroit et plus stable selon +les ages chez les non-fumeuses. Cela est cohérent avec la forme +sigmoïdale déjà repérée sur le diagramme en violon. +
+Cette analyse est à prendre avec des pincettes, c'est surtout pour +avoir un premier aperçu avant d'aller étudier et affiner les modèles +entraînés. C'est hélas la seule chose que nous permet la librairie +Seaborn puisque l'accès aux paramètres +ne +sera jamais permis! +
+ ++Construisons alors un modèle de classification binaire pour +chaque sous-groupe fumeuse/non fumeuse afin d'essayer d'évaluer +l'impact de la cigarette sur la capacité à avoir survécu pendant les +20 ans séparant les deux études: +
+from sklearn.model_selection import cross_val_score, cross_validate +from sklearn.model_selection import RepeatedStratifiedKFold + +def train(model, model_name, data): + """Entraine un modèle pour chaque partition(Fumeuse/Non Fumeuse). + Ajoute également 2 colonnes à data contenant les prédictions du modèles ainsi que les probabilités""" + + # On stocke les classifieurs pour chaque groupe (fumeuse/non fumeuse) + clf = dict(Yes=None, No=None, name=model_name) + + for is_smoker in ['Yes', 'No']: + # Définition des jeux de données pour la validation croisée + cv = RepeatedStratifiedKFold(n_splits=6, n_repeats=3, random_state=1) + X = data[data.Smoker == is_smoker]['Age'].to_frame() + targets = data[data.Smoker == is_smoker]['Death'] + + # Entrainement + scores = cross_validate(model, X, targets, scoring='roc_auc', cv=cv, n_jobs=-1, return_estimator=True) + # On récupère le meilleur classifieur + clf.update({is_smoker: scores['estimator'][scores['test_score'].argmax()]}) + + # Ajout des prédictions et probabilités du modèle entrainé pour la partition en cours(fumeuse ou non fumeuse) + # Les probabilités retournées sont [P(Death = 0), P(Death = 1)]. On stocke donc le deuxième élément. + # Visualiser ordre des classes{0,1} dans clf[is_smoker].classes_ pour retrouver cette information + data.loc[data.Smoker == is_smoker, f'DeathProba{model_name}'] = clf[is_smoker].predict_proba(data[data.Smoker == is_smoker]['Age'].to_frame())[:,1] + return clf ++
+Nous allons également entraîner un modèle d'arbre de décision, qui +peut se retrouver assez performant pour la classification binaire. +
+ +from sklearn.linear_model import LogisticRegression +from sklearn.tree import DecisionTreeClassifier + +# train models +model = LogisticRegression(solver='lbfgs', class_weight='balanced', C=0.1) +clf_LR = train(model, "LR", data) + +model = DecisionTreeClassifier(class_weight='balanced', min_samples_leaf=5) +clf_tree = train(model, "Tree", data) ++
+Observons les résultats +
+ +import numpy as np +from scipy.special import expit + +plt.clf() +colors = dict(Yes='green', No='Red') + +for is_smoker in ['Yes', 'No']: + clf = clf_LR[is_smoker] + age, death = data[data.Smoker == is_smoker][['Age', 'Death']].T.values + + # Prepare for drawing logistic curve estimate + age_values = np.linspace(-0, 100, 300) + logit = expit(age_values * clf.coef_ + clf.intercept_).ravel() + + # Draw + plt.scatter(age, death, color=colors[is_smoker], marker='+') + plt.plot(age_values, logit, color=colors[is_smoker], linewidth=2) + +figure_path = 'smoker_survival_logistic_custom.png' +plt.savefig(figure_path) +figure_path ++
+
+On a le même type de courbe que précèdemment. C'est encore
+contre-intuitif puisque les non-fumeuses(en vert) semblent se
+précipiter plus vite que les fumeuses (en rouge) vers la mort (passer de 0 à 1).
+
+Regardons quelques métriques usuelles +
+ +from sklearn import metrics + +metrics_dict = dict(accuracy=metrics.accuracy_score, + precision=metrics.precision_score, + recall=metrics.recall_score) + +for clf in [clf_LR, clf_tree]: + print(f"<<<clf is {clf['Yes'].__class__.__name__}>>>") + for is_smoker in ['Yes', 'No']: + print(f"------Smoker is {is_smoker}------") + y_true = data[data.Smoker == is_smoker]['Death'] + y_pred = clf[is_smoker].predict(data[data.Smoker == is_smoker]['Age'].to_frame()) + for metric,score in metrics_dict.items(): + print(f"{metric:>10} : {score(y_true, y_pred):.3}") + #print(metrics.confusion_matrix(y_true, y_pred))) + ++
+<<<clf is LogisticRegression>>> +------Smoker is Yes------ + accuracy : 0.722 + precision : 0.448 + recall : 0.712 +------Smoker is No------ + accuracy : 0.832 + precision : 0.683 + recall : 0.87 +<<<clf is DecisionTreeClassifier>>> +------Smoker is Yes------ + accuracy : 0.789 + precision : 0.537 + recall : 0.827 +------Smoker is No------ + accuracy : 0.851 + precision : 0.711 + recall : 0.887 ++ + +
from itertools import product +plt.clf() + +for clf,is_smoker in product([clf_LR, clf_tree], ['Yes', 'No']): + y_true = data[data.Smoker == is_smoker]['Death'] + y_pred_proba = data[data.Smoker == is_smoker][f"DeathProba{clf['name']}"] + + fpr, tpr, _ = metrics.roc_curve(y_true, y_pred_proba) + auc = metrics.roc_auc_score(y_true, y_pred_proba) + plt.plot(fpr,tpr,label=f"{clf['name']}(Smoker={is_smoker}), auc={auc:.2}") + plt.legend(loc=4) + +plt.xlabel('Taux de faux positifs') +plt.ylabel('Taux de vrais positifs') +plt.title('Courbe ROC') +figure_path = "roc.png" +plt.savefig(figure_path, bbox_inches='tight'); +plt.close() +return figure_path ++
+
+On se rend compte qu'il y a une différence non négligeable de +performance entre les modèles des fumeuses et ceux des non +fumeuses. Cette différence est visible même en changeant la nature du +modèle(Arbre de décision contre Régression logistique) +Avant de tirer des conclusions sur la nocivité du tabagisme, il +faudrait avoir des modèles ayant des performances similaires qui +disposeraient donc de capacité prédictives similaires. +
+ ++En conclusion, il faudrait plus de données pour avoir des +distributions plus représentatives de la réalité. Aussi, ajouter de +l'expertise métier(information a priori) aux modèles pourrait aider +avoir une bonne capacité prédictive et sortir du Paradox de Simpson. On imagine créer d'autres +variables, comme l'écart à l'espérance de vie par exemple. +
+