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

     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)

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.

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.