#+TITLE: Etude de la nocivité du tabac chez les femmes #+AUTHOR: Victor #+EMAIL:vtrltz@e.email * Préface :noexport: ** Propriétés du document #+LANGUAGE: French #+PROPERTY: header-args :session :exports both #+EXPORT_SELECT_TAGS: export #+EXPORT_EXCLUDE_TAGS: noexport #+bibliography: references.bib #+cite_export: basic author author-year ** Paramètres d'export #+HTML_HEAD: #+HTML_HEAD: #+HTML_HEAD: #+HTML_HEAD: #+HTML_HEAD: #+HTML_HEAD: ** Prérequis à l'utilisation de ce document Afin d'exécuter le document computationnel suivant, il convient de disposer des logiciels suivants : *** Emacs avec org-mode Afin d'exécuter les fragments de codes de ce document computationnel il convient d'utiliser une version d'Emacs intégrant Org-mode. *** Python = 3 #+BEGIN_SRC python :results output :exports both from sys import version_info if version_info.major < 3: print("Veuillez utiliser Python 3.") #+END_SRC #+RESULTS: #+BEGIN_SRC emacs-lisp :results output :exports both (unless (featurep 'ob-python) (print "Veuillez activer python dans org-babel (org-babel-do-languages) !")) #+END_SRC #+RESULTS: * Préambule ** Contexte de l'exercice Dans le cadre du [[https://www.fun-mooc.fr/fr/cours/recherche-reproductible-principes-methodologiques-pour-une-science-transparente/][/MOOC/ de l'INRIA sur la recherche reproductible]], je réalise un travail pratique évalué par les pairs. Ce travail pratique intervient au cours du module 3. Le sujet choisi est le numéro 6, intitulé : /Autour du Paradoxe de Simpson/. ** Origine des données Nous utilisons les données mises à disposition pour l'exercice par l'équipe pédagogique du MOOC. Ces données sont un extrait de celles d'études de l'incidence des pathologies de la thyroïde au sein de la population britannique. Furent exclus du jeu de donnée complet : - les hommes (n = 1285) ; - les femmes ayant arrêté de fumer (n = 162) - les femmes dont les données n'étaient pas disponibles (n = 18). URL des données : #+NAME: data-url https://gitlab.inria.fr/learninglab/mooc-rr/mooc-rr-ressources/-/raw/master/module3/Practical_session/Subject6_smoking.csv?inline=false Nous téléchargeons le jeu de donnée de l'exercice au format =.csv=. Dans ce jeu de données, chaque ligne correspond aux données d'une participante lors des études originales. Nous disposons des informations suivantes : | Nom de colonne | Description de la variable dépendante | |----------------+--------------------------------------------------------------------------------| | Smoker | Si la personne fume ou non. | | Status | Si la personne est vivante ou décédée. | | Age | L'âge de la personne si elle est vivante, l'âge à sa mort si elle est décédée. | ** Bibliothèques et fonctions utilisées Pour la suite des manipulations de données, nous faisons l'usage de plusieurs fonctions disponibles dans des bibliothèques Python. #+begin_src python :results output :exports code from os.path import exists from pandas import read_csv from matplotlib import pyplot as plt, rc from seaborn import histplot, kdeplot, lmplot from urllib import request from numbers import Real from numpy import nan, arange #+end_src #+RESULTS: ** Paramètres communs des graphiques Nous souhaitons que les graphiques soient réalisés avec TeX. #+begin_src python :results output :exports code # Draw in LaTeX rc('font', **{'family': 'serif', 'serif': ['Courier New']}) plt.rcParams['text.usetex'] = True #+end_src #+RESULTS: ** Import des données Dans le cas où une copie locale n'est pas présente, nous téléchargeons le jeu de données complet depuis son adresse de dépôt. #+begin_src python :results output :exports both :var data_url=data-url data_file = "Subject6_smoking.csv" if not exists(data_file): request.urlretrieve(data_url, data_file) print('Data set downloaded. \n') data = read_csv(filepath_or_buffer='Subject6_smoking.csv') data.info(verbose=True) print('\n', data.head()) #+end_src #+RESULTS: #+begin_example RangeIndex: 1314 entries, 0 to 1313 Data columns (total 3 columns): Smoker 1314 non-null object Status 1314 non-null object Age 1314 non-null float64 dtypes: float64(1), object(2) memory usage: 30.9+ KB 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 #+end_example ** Vérifications automatiques des données Vérifions le jeu de données : - nous ne devons avoir que les valeurs /Yes/ et /No/ dans la colonne /Smoker/ ; - seules les catégories /ALive/ et /Dead/ doivent être présentent dans la colonne /Status/ ; - l'âge doit être un réel entre 0 et 123 (le record de longévité attesté étant de moins de 123 ans). En cas de valeur incorrecte, nous remplaçons ladite valeur par =nan= (/Not a number/). #+begin_src python :results output :exports code for row in data.iterrows(): if row[1][0] != 'Yes' and row[1][0] != 'No': print(f"Uncorrect value in 'Smoker' column: {row[1][0]} (index = {row[0]})") data['Smoker'].iat[row[0]] = nan if row[1][1] != 'Alive' and row[1][1] != 'Dead': print(f"Uncorrect value in 'Status' column: {row[1][1]} (index = {row[0]})") data['Status'].iat[row[0]] = nan if not isinstance(row[1][2], Real) and not 0 <= row[1][2] <= 123: print(f"Uncorrect value in 'Age' column: {row[1][2]} (index = {row[0]})") data['Age'].iat[row[0]] = nan #+end_src #+RESULTS: ** Visualisations des données brutes *** Distribution des âges des participantes Nous souhaitons visualiser les distributions des âges selon les habitudes de tabagisme. #+begin_src python :results output file :exports none fig, ax = plt.subplots() kdeplot(data=data, x=data['Age'], hue=data['Smoker']) plt.savefig("group_ages_kde.png") plt.close() #+end_src #+RESULTS: [[file:]] #+CAPTION: Distribution des âges des participantes selon les habitudes de tabagisme (estimation par noyau) #+NAME: fig:age-kde [[./group_ages_kde.png]] Constatons que les distributions en âge des participantes selon les habitudes de tabagisme ne correspondent pas : les femmes de plus de 60 ans sont davantage présentes dans le groupe des non-fumeuses. * Analyses de la nocivité du tabagisme chez les femmes Dans les analyses suivantes, nous souhaitons évaluer la nocivité du tabagisme chez les femmes. Dans le cadre de notre exercice, nous utilisons trois approches différentes. Notre hypothèse est que le tabac est nocif chez les femmes. ** Approche 1 | Taux de mortalité des fumeuses et non-fumeuses Nous cherchons à déterminer les taux de mortalité des deux groupes de notre étude. Nous souhaitons associer ces taux de mortalité à des intervalles de confiance à 95%. *** Calcul des paramètres Commençons par séparer les deux groupes. Le nombre de participantes mortes dans un groupe ramené à l'effectif total du groupe nous donne le taux de mortalité dudit groupe. #+begin_src python :results output :exports code rounding = 1 smokers = data[(data['Smoker'] == 'Yes')] nonsmokers = data[(data['Smoker'] == 'No')] dead_smokers = smokers[(smokers['Status'] == 'Dead')] dead_nonsmokers = nonsmokers[(nonsmokers['Status'] == 'Dead')] smokers_death_rate = round(dead_smokers.shape[0] / smokers.shape[0] * 100, rounding) nonsmokers_death_rate = round(dead_nonsmokers.shape[0] / nonsmokers.shape[0] * 100, rounding) #+end_src #+RESULTS: Nous déterminons l'intervalle de confiance à 95% avec la formule suivante : #+begin_latex $\text{IC}_{95} = \frac{100}{n} d \pm 1.96 \sqrt{d}$ #+end_latex Où =n= correspond à l'effectif du groupe considéré et =d= au nombre d'individus décédés au sein du même groupe. Voici le code que nous avons utilisé pour réaliser le calcul des intervalles de confiance : #+begin_src python :results output :exports code def death_rate_ci95(sample_size: int, deads: int, per_people: int = 100) -> tuple: """Compute the 95% confidence interval upper and lower limits of a mortality rate in a sample. :param sample_size: number of people in the sample. :param deads: number of deads people in the sample. :param per_people: standardisation number: deads for per_people people. """ from math import sqrt upper = per_people / sample_size * (deads + 1.96 * sqrt(deads)) lower = per_people / sample_size * (deads - 1.96 * sqrt(deads)) return lower, upper smokers_dr_95ic = death_rate_ci95(smokers.shape[0], dead_smokers.shape[0]) nonsmokers_dr_95ic = death_rate_ci95(nonsmokers.shape[0], dead_nonsmokers.shape[0]) #+end_src #+RESULTS: *** Résultats #+begin_src python :results output table raw :exports results print(f"|Tabagisme|Effectif total|Individus décédés|Taux de mortalité (%)|IC95%|\n" f"|-+-+-+-+-|\n" f"|Fumeuses|{smokers.shape[0]}|{dead_smokers.shape[0]}|{smokers_death_rate}|[{round(smokers_dr_95ic[0], rounding)};{round(smokers_dr_95ic[1], rounding)}]|\n" f"|Non fumeuses|{nonsmokers.shape[0]}|{dead_nonsmokers.shape[0]}|{nonsmokers_death_rate}|[{round(nonsmokers_dr_95ic[0], rounding)};{round(nonsmokers_dr_95ic[1], rounding)}]|") #+end_src #+RESULTS: | Tabagisme | Effectif total | Individus décédés | Taux de mortalité (%) | IC95% | |--------------+----------------+-------------------+-----------------------+-------------| | Fumeuses | 582 | 139 | 23.9 | [19.9;27.9] | | Non fumeuses | 732 | 230 | 31.4 | [27.4;35.5] | #+begin_src python :results output file :exports none fig, ax = plt.subplots(nrows=1, ncols=2, sharex=True) histplot(data=data, x='Smoker', hue='Status', multiple='stack', shrink=0.8, ax=ax[0]) ax[1].bar(x=[0, 1], height=[smokers_death_rate, nonsmokers_death_rate], yerr=[smokers_dr_95ic[1] - smokers_death_rate, nonsmokers_dr_95ic[1] - nonsmokers_death_rate], color='grey', edgecolor='k', error_kw=dict(lolims=True)) # In order to modify the default lolims error bars shape for ch in ax[1].get_children(): if str(ch).startswith('Line2D'): ch.set_marker('_') ch.set_markersize(6) ax[1].set_ylabel(r'Death rate (\%)') ax[1].set_xlabel('Smoker') plt.tight_layout() plt.savefig("groups_death_rate.png") #+end_src #+RESULTS: [[file:]] #+CAPTION: Effectifs et taux de mortalité parmi les fumeuses et les non fumeuses. #+NAME: fig:effectifs-mortalite [[./groups_death_rate.png]] Pour le jeu de donnée considéré, le taux de mortalité des fumeuses est de src_python[:exports results]{smokers_death_rate} {{{results(=23.88=)}}} contre src_python[:exports results]{nonsmokers_death_rate} {{{results(=31.42=)}}} pour celui observé chez les non-fumeuses. Ce qui va à l'encontre de notre hypothèse initiale. Il est possible que l'hétérogénéité constatée dans la distribution des âges explique ce résultat contre intuitif. ** Approche 2 | Taux de mortalité des fumeuses et non-fumeuses selon la classe d'âge Nous étudions désormais le taux de mortalité selon les habitudes de tabagisme et selon la classe d'âge. Nous considérons quatre classes : - [18 ; 34[ - [34 ; 54[ - [54 ; 64[ - [64 ; 123[ *** Calcul des paramètres et résultats Premièrement nous devons séparer le jeu de donnée selon les classes d'âge. Puis nous procédons de la même manière que précédemment : nous identifions les effectifs de chaque groupe, le nombre de décès au sein des sous-groupes et nous calculons le taux de mortalité associé à son intervalle de confiance à 95%. #+begin_src python :results output table raw :exports both categories = [(18, 34), (34, 54), (54, 64), (64, 123)] death_rates = [] # plotting purpose dr_uncertainties = [] # plotting purpose print(f"|Classe d'âge|Tabagisme|Effectif total|Individus décédés|\ Taux de mortalité (%)|IC95%|\n" f"|-+-+-+-+-+-|") for category in categories: for smoker in data['Smoker'].unique(): alive, dead = nan, nan if smoker == 'Yes': smoke = 'Fumeuses' else: smoke = 'Non-fumeuses' for status in data['Status'].unique(): subset = data[(data['Age'] >= category[0]) & (data['Age'] < category[1]) & (data['Smoker'] == smoker) & (data['Status'] == status)] if status == 'Alive': alive = subset.shape[0] else: dead = subset.shape[0] death_rate = dead / (alive + dead) * 100 death_rates.append(death_rate) dr_ci95 = death_rate_ci95(alive + dead, dead) dr_uncertainties.append(dr_ci95[1] - death_rate) rounded_dr_ci95 = list(dr_ci95) rounded_dr_ci95[0] = round(rounded_dr_ci95[0], rounding) rounded_dr_ci95[1] = round(rounded_dr_ci95[1], rounding) r_dr_ci95 = tuple(rounded_dr_ci95) print(f"|[{category[0]};{category[1]}[|{smoke}|{alive + dead}|\ {dead}|{round(death_rate, rounding)}|{r_dr_ci95}|") #+end_src #+RESULTS: | Classe d'âge | Tabagisme | Effectifs total | Individus décédés | Taux de mortalité (%) | IC95% | |--------------+--------------+-----------------+-------------------+-----------------------+---------------| | [18;34[ | Fumeuses | 179 | 5 | 2.8 | (0.3, 5.2) | | [18;34[ | Non fumeuses | 219 | 6 | 2.7 | (0.5, 4.9) | | [34;54[ | Fumeuses | 239 | 41 | 17.2 | (11.9, 22.4) | | [34;54[ | Non fumeuses | 199 | 19 | 9.5 | (5.3, 13.8) | | [54;64[ | Fumeuses | 115 | 51 | 44.3 | (32.2, 56.5) | | [54;64[ | Non fumeuses | 119 | 39 | 32.8 | (22.5, 43.1) | | [64;123[ | Fumeuses | 49 | 42 | 85.7 | (59.8, 111.6) | | [64;123[ | Non fumeuses | 195 | 166 | 85.1 | (72.2, 98.1) | #+begin_src python :results output file :exports both fig, ax = plt.subplots() x_pos = arange(len(categories)) width = 0.35 xlabels = list(map(str, categories)) for x_l, xlabel in enumerate(xlabels): head, _, tail = xlabel.partition(',') xlabels[x_l] = '[' + head[1:] + ';' + tail[:-1] + '[' plt.bar(x=x_pos - width / 2, height=death_rates[0:len(death_rates):2], yerr=dr_uncertainties[0:len(death_rates):2], label='Smoker', width=width, edgecolor='k', error_kw=dict(lolims=True)) plt.bar(x=x_pos + width / 2, height=death_rates[1:len(death_rates):2], yerr=dr_uncertainties[1:len(death_rates):2], label='Nonsmoker', width=width, edgecolor='k', error_kw=dict(lolims=True)) # In order to modify the default lolims error bars shape for ch in ax.get_children(): if str(ch).startswith('Line2D'): ch.set_marker('_') ch.set_markersize(6) plt.ylabel('Death rate (%)') plt.xlabel('Age category') plt.xticks(x_pos) ax.set_xticklabels(xlabels) plt.legend() plt.tight_layout() plt.savefig("age_categories_death_rate.png") #+end_src #+RESULTS: [[file:]] #+CAPTION: Taux de mortalité et IC95% selon les habitudes de tabagisme et la classe d'âge #+NAME: fig:mortalite-classe-age [[./age_categories_death_rate.png]] Pour chacune de nos classes d'âge, les résultats indiquent que le taux de mortalité des fumeuses est supérieur à celui des non-fumeuses ; soit une interprétation des données en contradiction avec la première analyse sans prise en compte de l'âge. Outre un biais potentiel lié aux distributions des âges hétérogènes, d'autres biais sont susceptibles d'exister dans les analyses précédentes : - un biais lié à la non prise en compte de l'âge dans le calcul des taux de mortalité ; - un biais lié aux choix de catégorisation. En somme au moins trois biais peuvent se combiner dans nos études précédentes et fausser nos interprétations. Enfin, notons que l'âge semble avoir un effet supérieur au tabac sur le taux de mortalité puisque l'amplitude des différences entre les classes d'âge est supérieure à celle entre les habitudes de tabagisme. En somme : le tabac tue, mais vieillir tue davantage. ** Approche 3 | Probabilité de décès des fumeuses et non-fumeuses selon l'âge Evitons les deux derniers biais en passant du calcul du taux de mortalité au calcul de la probabilité de décès selon l'âge et les habitudes de tabagisme. Pour ce faire nous pouvons utiliser la régression logistique. *** Manipulation des données Plutôt que de disposer des status =Alive= et =Dead=, nous voulons qu'ils soient respectivement codés avec =0= et =1=. Nous souhaitons à nouveau disposer de jeux de données spécifiques à chaque groupe : fumeuses et non fumeuses. #+begin_src python :results output :exports both deaths = [] for status in data['Status']: if status == 'Alive': deaths.append(0) else: deaths.append(1) data['Death'] = deaths smokers = data[(data['Smoker'] == 'Yes')] nonsmokers = data[(data['Smoker'] == 'No')] print(smokers.head(), '\n', '\n', nonsmokers.head()) #+end_src #+RESULTS: #+begin_example Smoker Status Age Death 0 Yes Alive 21.0 0 1 Yes Alive 19.3 0 4 Yes Alive 81.4 0 7 Yes Dead 57.5 1 8 Yes Alive 24.8 0 Smoker Status Age Death 2 No Dead 57.5 1 3 No Alive 47.1 0 5 No Alive 36.8 0 6 No Alive 23.8 0 11 No Dead 66.0 1 #+end_example *** Régressions logistiques **** Réglages et code de la modélisation Pour effectuer la régression logistique, nous utilisons la bibliothèque =statsmodels=. Nous laissons les paramètres par défaut : le modèle est le logit, la méthode est celle du maximum de vraisemblance (/Maximum Likelyhood Estimation/). #+begin_src python :results output :exports code def logistic_regression(independent_variable: list, dependent_variable: list) -> object: """Logisitic regression with Logit model and MLE method. :param independent_variable: the independent variable :param dependent_variable: the dependent variable :return fitted_model: the logisitc regression fitted model. """ import statsmodels.api as sm model = sm.Logit(dependent_variable, sm.add_constant(independent_variable)) fitted_model = model.fit() print('\n', fitted_model.summary()) return fitted_model #+end_src #+RESULTS: **** Cas des fumeuses #+begin_src python :results output :exports both proba_smokers, model_smokers = logistic_regression(smokers['Age'], smokers['Death']) #+end_src #+RESULTS: #+begin_example /usr/lib/python3/dist-packages/numpy/core/fromnumeric.py:2495: FutureWarning: Method .ptp is deprecated and will be removed in a future version. Use numpy.ptp instead. return ptp(axis=axis, out=out, **kwargs) Optimization terminated successfully. Current function value: 0.412727 Iterations 7 Logit Regression Results ============================================================================== Dep. Variable: Death No. Observations: 582 Model: Logit Df Residuals: 580 Method: MLE Df Model: 1 Date: jeu., 14 oct. 2021 Pseudo R-squ.: 0.2492 Time: 05:32:17 Log-Likelihood: -240.21 converged: True LL-Null: -319.94 Covariance Type: nonrobust LLR p-value: 1.477e-36 ============================================================================== coef std err z P>|z| [0.025 0.975] ------------------------------------------------------------------------------ const -5.5081 0.466 -11.814 0.000 -6.422 -4.594 Age 0.0890 0.009 10.203 0.000 0.072 0.106 ============================================================================== #+end_example **** Cas des non-fumeuses #+begin_src python :results output :exports both proba_nonsmokers, model_nonsmokers = logistic_regression(nonsmokers['Age'], nonsmokers['Death']) #+end_src #+RESULTS: #+begin_example /usr/lib/python3/dist-packages/numpy/core/fromnumeric.py:2495: FutureWarning: Method .ptp is deprecated and will be removed in a future version. Use numpy.ptp instead. return ptp(axis=axis, out=out, **kwargs) Optimization terminated successfully. Current function value: 0.354560 Iterations 7 Logit Regression Results ============================================================================== Dep. Variable: Death No. Observations: 732 Model: Logit Df Residuals: 730 Method: MLE Df Model: 1 Date: jeu., 14 oct. 2021 Pseudo R-squ.: 0.4304 Time: 05:32:21 Log-Likelihood: -259.54 converged: True LL-Null: -455.62 Covariance Type: nonrobust LLR p-value: 2.808e-87 ============================================================================== coef std err z P>|z| [0.025 0.975] ------------------------------------------------------------------------------ const -6.7955 0.479 -14.174 0.000 -7.735 -5.856 Age 0.1073 0.008 13.742 0.000 0.092 0.123 ============================================================================== #+end_example *** Visualisation des résultats #+begin_src python :results output file :exports both fig, ax = plt.subplots() ax.set_xlim(15, 100) ax.set_ylabel('Death probability') ax.set_xlabel('Age') lmplot(data=data, x='Age', y='Death', hue='Smoker', logistic=True, y_jitter=0.01, truncate=False, markers=['o', 'x'], scatter_kws=dict(alpha=0.3)) plt.tight_layout() plt.savefig('death_probabilities.png') #+end_src #+RESULTS: [[file:]] #+CAPTION: Estimation des probabilités de décès selon l'âge et les habitudes de tabagisme. #+NAME: fig:proba-deces [[./death_probabilities.png]] L'analyse graphique des résultats semble indiquer une plus forte probabilité de décès des femmes fumeuses avant 55 ans. Après cet âge, les probabilités de décès ne paraissent plus influencées par les habitudes de tabagisme. Le tabac chez les femmes est donc nocif, surtout avant 55 ans. * Paradoxe de Simpson Nous sommes face à un paradoxe de Simpson quand la tendance observée au sein de sous-groupes disparaît ou est inversée quand les sous-groupes sont combinés. ** Qu'en est-il du paradoxe de Simpson dans notre cas ? Lors de notre première approche, les résultats indiquent que le taux de mortalité est supérieur chez les non-fumeuses (src_python[:exports results]{nonsmokers_death_rate} {{{results(=31.42=)}}} contre src_python[:exports results]{smokers_death_rate}{{{results(=23.88=)}}}). Lors de la deuxième approche, lors de la prise en compte de l'âge, les résultats indiquent l'inverse : le taux de mortalité est supérieur chez les fumeuses. Nous sommes face à un paradoxe de Simpson : la prise en compte ou non de la variable =Age= inverse les interprétations que nous pouvons faire des effets du tabagisme chez les femmes. ** Quelles sont les explications dans notre cas ? Il est probable que le paradoxe observé dans le cadre de cet exercice s'explique en partie par les faits suivants : - il existe une différence de src_python[:exports results]{nonsmokers.shape[0] - smokers.shape[0]} {{{results(=150=)}}} individus entre les effectifs des fumeuses et des non fumeuses, soit environ src_python[:exports results]{round((nonsmokers.shape[0] - smokers.shape[0])/(nonsmokers.shape[0] + smokers.shape[0]) * 100, 0)} {{{results(=11.0=)}}} % de l'effectif total, ce qui n'est pas négligeable ; - les distributions en âge au sein des groupes (selon les habitudes de tabagisme) ne coïncident pas ; - l'âge a un effet supérieur au tabac sur le taux de mortalité.