Sujet 6: Autour du Paradoxe de Simpson

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

data.describe(include='all')

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)

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

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

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

[[file:[2. 1. 0.92939002 0.91644794 0.89769684 0.8143829 0.76693106 0.7327737 0.68688525 0.63702014 0.59390503 0.52309613 0.46737312 0.42237903 0.35418428 0.21520288 0. ]]]

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)

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

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)

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

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

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

[[file:[2. 1. 0.92939002 0.91644794 0.89769684 0.8143829 0.76693106 0.7327737 0.68688525 0.63702014 0.59390503 0.52309613 0.46737312 0.42237903 0.35418428 0.21520288 0. ]]] 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.

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

[[file:[2. 1. 0.92939002 0.91644794 0.89769684 0.8143829 0.76693106 0.7327737 0.68688525 0.63702014 0.59390503 0.52309613 0.46737312 0.42237903 0.35418428 0.21520288 0. ]]]

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']]

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

[[file:[2. 1. 0.92939002 0.91644794 0.89769684 0.8143829 0.76693106 0.7327737 0.68688525 0.63702014 0.59390503 0.52309613 0.46737312 0.42237903 0.35418428 0.21520288 0. ]]]

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

[[file:[2. 1. 0.92939002 0.91644794 0.89769684 0.8143829 0.76693106 0.7327737 0.68688525 0.63702014 0.59390503 0.52309613 0.46737312 0.42237903 0.35418428 0.21520288 0. ]]] 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

[[file:[2. 1. 0.92939002 0.91644794 0.89769684 0.8143829 0.76693106 0.7327737 0.68688525 0.63702014 0.59390503 0.52309613 0.46737312 0.42237903 0.35418428 0.21520288 0. ]]] [[file:[2. 1. 0.92939002 0.91644794 0.89769684 0.8143829 0.76693106 0.7327737 0.68688525 0.63702014 0.59390503 0.52309613 0.46737312 0.42237903 0.35418428 0.21520288 0. ]]]

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.