Commit c31a1e51 authored by x0s's avatar x0s

Ajout sujet 6 paradoxe simpson affichage okay mais en trichant

parent 5eb6ca4e
<?xml version="1.0" encoding="utf-8"?>
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN"
"http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
<html xmlns="http://www.w3.org/1999/xhtml" lang="fr" xml:lang="fr">
<head>
<!-- 2020-04-24 ven. 09:47 -->
<meta http-equiv="Content-Type" content="text/html;charset=utf-8" />
<meta name="viewport" content="width=device-width, initial-scale=1" />
<title>Sujet 6: Autour du Paradoxe de Simpson</title>
<meta name="generator" content="Org mode" />
<meta name="author" content="Guillaume" />
<style type="text/css">
<!--/*--><![CDATA[/*><!--*/
.title { text-align: center;
margin-bottom: .2em; }
.subtitle { text-align: center;
font-size: medium;
font-weight: bold;
margin-top:0; }
.todo { font-family: monospace; color: red; }
.done { font-family: monospace; color: green; }
.priority { font-family: monospace; color: orange; }
.tag { background-color: #eee; font-family: monospace;
padding: 2px; font-size: 80%; font-weight: normal; }
.timestamp { color: #bebebe; }
.timestamp-kwd { color: #5f9ea0; }
.org-right { margin-left: auto; margin-right: 0px; text-align: right; }
.org-left { margin-left: 0px; margin-right: auto; text-align: left; }
.org-center { margin-left: auto; margin-right: auto; text-align: center; }
.underline { text-decoration: underline; }
#postamble p, #preamble p { font-size: 90%; margin: .2em; }
p.verse { margin-left: 3%; }
pre {
border: 1px solid #ccc;
box-shadow: 3px 3px 3px #eee;
padding: 8pt;
font-family: monospace;
overflow: auto;
margin: 1.2em;
}
pre.src {
position: relative;
overflow: visible;
padding-top: 1.2em;
}
pre.src:before {
display: none;
position: absolute;
background-color: white;
top: -10px;
right: 10px;
padding: 3px;
border: 1px solid black;
}
pre.src:hover:before { display: inline;}
/* Languages per Org manual */
pre.src-asymptote:before { content: 'Asymptote'; }
pre.src-awk:before { content: 'Awk'; }
pre.src-C:before { content: 'C'; }
/* pre.src-C++ doesn't work in CSS */
pre.src-clojure:before { content: 'Clojure'; }
pre.src-css:before { content: 'CSS'; }
pre.src-D:before { content: 'D'; }
pre.src-ditaa:before { content: 'ditaa'; }
pre.src-dot:before { content: 'Graphviz'; }
pre.src-calc:before { content: 'Emacs Calc'; }
pre.src-emacs-lisp:before { content: 'Emacs Lisp'; }
pre.src-fortran:before { content: 'Fortran'; }
pre.src-gnuplot:before { content: 'gnuplot'; }
pre.src-haskell:before { content: 'Haskell'; }
pre.src-hledger:before { content: 'hledger'; }
pre.src-java:before { content: 'Java'; }
pre.src-js:before { content: 'Javascript'; }
pre.src-latex:before { content: 'LaTeX'; }
pre.src-ledger:before { content: 'Ledger'; }
pre.src-lisp:before { content: 'Lisp'; }
pre.src-lilypond:before { content: 'Lilypond'; }
pre.src-lua:before { content: 'Lua'; }
pre.src-matlab:before { content: 'MATLAB'; }
pre.src-mscgen:before { content: 'Mscgen'; }
pre.src-ocaml:before { content: 'Objective Caml'; }
pre.src-octave:before { content: 'Octave'; }
pre.src-org:before { content: 'Org mode'; }
pre.src-oz:before { content: 'OZ'; }
pre.src-plantuml:before { content: 'Plantuml'; }
pre.src-processing:before { content: 'Processing.js'; }
pre.src-python:before { content: 'Python'; }
pre.src-R:before { content: 'R'; }
pre.src-ruby:before { content: 'Ruby'; }
pre.src-sass:before { content: 'Sass'; }
pre.src-scheme:before { content: 'Scheme'; }
pre.src-screen:before { content: 'Gnu Screen'; }
pre.src-sed:before { content: 'Sed'; }
pre.src-sh:before { content: 'shell'; }
pre.src-sql:before { content: 'SQL'; }
pre.src-sqlite:before { content: 'SQLite'; }
/* additional languages in org.el's org-babel-load-languages alist */
pre.src-forth:before { content: 'Forth'; }
pre.src-io:before { content: 'IO'; }
pre.src-J:before { content: 'J'; }
pre.src-makefile:before { content: 'Makefile'; }
pre.src-maxima:before { content: 'Maxima'; }
pre.src-perl:before { content: 'Perl'; }
pre.src-picolisp:before { content: 'Pico Lisp'; }
pre.src-scala:before { content: 'Scala'; }
pre.src-shell:before { content: 'Shell Script'; }
pre.src-ebnf2ps:before { content: 'ebfn2ps'; }
/* additional language identifiers per "defun org-babel-execute"
in ob-*.el */
pre.src-cpp:before { content: 'C++'; }
pre.src-abc:before { content: 'ABC'; }
pre.src-coq:before { content: 'Coq'; }
pre.src-groovy:before { content: 'Groovy'; }
/* additional language identifiers from org-babel-shell-names in
ob-shell.el: ob-shell is the only babel language using a lambda to put
the execution function name together. */
pre.src-bash:before { content: 'bash'; }
pre.src-csh:before { content: 'csh'; }
pre.src-ash:before { content: 'ash'; }
pre.src-dash:before { content: 'dash'; }
pre.src-ksh:before { content: 'ksh'; }
pre.src-mksh:before { content: 'mksh'; }
pre.src-posh:before { content: 'posh'; }
/* Additional Emacs modes also supported by the LaTeX listings package */
pre.src-ada:before { content: 'Ada'; }
pre.src-asm:before { content: 'Assembler'; }
pre.src-caml:before { content: 'Caml'; }
pre.src-delphi:before { content: 'Delphi'; }
pre.src-html:before { content: 'HTML'; }
pre.src-idl:before { content: 'IDL'; }
pre.src-mercury:before { content: 'Mercury'; }
pre.src-metapost:before { content: 'MetaPost'; }
pre.src-modula-2:before { content: 'Modula-2'; }
pre.src-pascal:before { content: 'Pascal'; }
pre.src-ps:before { content: 'PostScript'; }
pre.src-prolog:before { content: 'Prolog'; }
pre.src-simula:before { content: 'Simula'; }
pre.src-tcl:before { content: 'tcl'; }
pre.src-tex:before { content: 'TeX'; }
pre.src-plain-tex:before { content: 'Plain TeX'; }
pre.src-verilog:before { content: 'Verilog'; }
pre.src-vhdl:before { content: 'VHDL'; }
pre.src-xml:before { content: 'XML'; }
pre.src-nxml:before { content: 'XML'; }
/* add a generic configuration mode; LaTeX export needs an additional
(add-to-list 'org-latex-listings-langs '(conf " ")) in .emacs */
pre.src-conf:before { content: 'Configuration File'; }
table { border-collapse:collapse; }
caption.t-above { caption-side: top; }
caption.t-bottom { caption-side: bottom; }
td, th { vertical-align:top; }
th.org-right { text-align: center; }
th.org-left { text-align: center; }
th.org-center { text-align: center; }
td.org-right { text-align: right; }
td.org-left { text-align: left; }
td.org-center { text-align: center; }
dt { font-weight: bold; }
.footpara { display: inline; }
.footdef { margin-bottom: 1em; }
.figure { padding: 1em; }
.figure p { text-align: center; }
.inlinetask {
padding: 10px;
border: 2px solid gray;
margin: 10px;
background: #ffffcc;
}
#org-div-home-and-up
{ text-align: right; font-size: 70%; white-space: nowrap; }
textarea { overflow-x: auto; }
.linenr { font-size: smaller }
.code-highlighted { background-color: #ffff00; }
.org-info-js_info-navigation { border-style: none; }
#org-info-js_console-label
{ font-size: 10px; font-weight: bold; white-space: nowrap; }
.org-info-js_search-highlight
{ background-color: #ffff00; color: #000000; font-weight: bold; }
.org-svg { width: 90%; }
/*]]>*/-->
</style>
<link rel="stylesheet" type="text/css" href="http://www.pirilampo.org/styles/readtheorg/css/htmlize.css"/>
<link rel="stylesheet" type="text/css" href="http://www.pirilampo.org/styles/readtheorg/css/readtheorg.css"/>
<script src="https://ajax.googleapis.com/ajax/libs/jquery/2.1.3/jquery.min.js"></script>
<script src="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.4/js/bootstrap.min.js"></script>
<script type="text/javascript" src="http://www.pirilampo.org/styles/lib/js/jquery.stickytableheaders.js"></script>
<script type="text/javascript" src="http://www.pirilampo.org/styles/readtheorg/js/readtheorg.js"></script>
<script type="text/javascript">
/*
@licstart The following is the entire license notice for the
JavaScript code in this tag.
Copyright (C) 2012-2019 Free Software Foundation, Inc.
The JavaScript code in this tag is free software: you can
redistribute it and/or modify it under the terms of the GNU
General Public License (GNU GPL) as published by the Free Software
Foundation, either version 3 of the License, or (at your option)
any later version. The code is distributed WITHOUT ANY WARRANTY;
without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU GPL for more details.
As additional permission under GNU GPL version 3 section 7, you
may distribute non-source (e.g., minimized or compacted) forms of
that code without the copy of the GNU GPL normally required by
section 4, provided you include this license notice and a URL
through which recipients can access the Corresponding Source.
@licend The above is the entire license notice
for the JavaScript code in this tag.
*/
<!--/*--><![CDATA[/*><!--*/
function CodeHighlightOn(elem, id)
{
var target = document.getElementById(id);
if(null != target) {
elem.cacheClassElem = elem.className;
elem.cacheClassTarget = target.className;
target.className = "code-highlighted";
elem.className = "code-highlighted";
}
}
function CodeHighlightOff(elem, id)
{
var target = document.getElementById(id);
if(elem.cacheClassElem)
elem.className = elem.cacheClassElem;
if(elem.cacheClassTarget)
target.className = elem.cacheClassTarget;
}
/*]]>*///-->
</script>
<script type="text/x-mathjax-config">
MathJax.Hub.Config({
displayAlign: "center",
displayIndent: "0em",
"HTML-CSS": { scale: 100,
linebreaks: { automatic: "false" },
webFont: "TeX"
},
SVG: {scale: 100,
linebreaks: { automatic: "false" },
font: "TeX"},
NativeMML: {scale: 100},
TeX: { equationNumbers: {autoNumber: "AMS"},
MultLineWidth: "85%",
TagSide: "right",
TagIndent: ".8em"
}
});
</script>
<script type="text/javascript"
src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.0/MathJax.js?config=TeX-AMS_HTML"></script>
</head>
<body>
<div id="content">
<h1 class="title">Sujet 6: Autour du Paradoxe de Simpson</h1>
<div id="table-of-contents">
<h2>Table des matières</h2>
<div id="text-table-of-contents">
<ul>
<li><a href="#org4c85e26">1. Introduction</a></li>
<li><a href="#orgf6f7892">2. Mission</a>
<ul>
<li><a href="#orge971908">2.1. Relation entre tabagisme et décès</a></li>
<li><a href="#org3788f96">2.2. Partitionnement selon la classe d'âge</a></li>
<li><a href="#org25b0d41">2.3. Modélisation de la probabilité de décès en fonction de l'âge</a></li>
</ul>
</li>
</ul>
</div>
</div>
<div id="outline-container-org4c85e26" class="outline-2">
<h2 id="org4c85e26"><span class="section-number-2">1</span> Introduction</h2>
<div class="outline-text-2" id="text-1">
<p>
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.
</p>
<p>
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.
(<a href="https://gitlab.inria.fr/learninglab/mooc-rr/mooc-rr-ressources/blob/master//module3/Practical_session/Subject6_fr.md">source</a>)
</p>
</div>
</div>
<div id="outline-container-orgf6f7892" class="outline-2">
<h2 id="orgf6f7892"><span class="section-number-2">2</span> Mission</h2>
<div class="outline-text-2" id="text-2">
</div>
<div id="outline-container-orge971908" class="outline-3">
<h3 id="orge971908"><span class="section-number-3">2.1</span> Relation entre tabagisme et décès</h3>
<div class="outline-text-3" id="text-2-1">
<p>
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 ?
</p>
<p>
Tout d'abord Chargons le jeu de données dans une dataframe pandas pour
faciliter la manipulation de données
</p>
<div class="org-src-container">
<pre class="src src-python"><span style="color: #a020f0;">from</span> pathlib <span style="color: #a020f0;">import</span> Path
<span style="color: #a020f0;">import</span> matplotlib.pyplot <span style="color: #a020f0;">as</span> plt
<span style="color: #a020f0;">import</span> pandas <span style="color: #a020f0;">as</span> pd
<span style="color: #b22222;"># </span><span style="color: #b22222;">set floating precision to 1 (ie: 97.244-&gt;97.2) </span>
pd.set_option(<span style="color: #8b2252;">"display.precision"</span>, 1)
<span style="color: #a0522d;">data_url</span> = <span style="color: #8b2252;">"http://gitlab.inria.fr/learninglab/mooc-rr/mooc-rr-ressources/blob/master/module3/Practical_session/Subject6_smoking.csv"</span>
<span style="color: #a0522d;">data_file</span> = Path(<span style="color: #8b2252;">"Subject6_smoking.csv"</span>)
<span style="color: #a020f0;">if</span> <span style="color: #a020f0;">not</span> data_file.exists():
<span style="color: #a020f0;">raise</span> <span style="color: #228b22;">IOError</span>(f<span style="color: #8b2252;">"Jeu de donn&#233;es introuvable, merci de le t&#233;l&#233;charger &#224; nouveau via ce lien {data_url}"</span>)
<span style="color: #a0522d;">data</span> = pd.read_csv(data_file, dtype={<span style="color: #8b2252;">'Smoker'</span>: <span style="color: #8b2252;">'category'</span>, <span style="color: #8b2252;">'Status'</span>: <span style="color: #8b2252;">'category'</span>})
data
</pre>
</div>
<pre class="example">
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]
</pre>
<div class="org-src-container">
<pre class="src src-python">data.describe(include=<span style="color: #8b2252;">'all'</span>)
</pre>
</div>
<pre class="example">
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
</pre>
<p>
On a bien deux valeurs possibles pour les variables catégorielles:
</p>
<ul class="org-ul">
<li>Fumeuse(Smoker): Yes/No
<ul class="org-ul">
<li>Indique si la personne est fumeuse <b>au moment de la seconde étude</b>
(20 ans après la première)</li>
</ul></li>
<li>Etat(Status) : Alive/Dead
<ul class="org-ul">
<li>Indique si la personne est vivante <b>au moment de la seconde étude</b></li>
</ul></li>
</ul>
<p>
<b>Au moment de la première étude</b>, 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.
</p>
<p>
Etudions le lien entre tabagisme et décès via ce tableau de
contingence sur les variables Smoker et Status:
</p>
<div class="org-src-container">
<pre class="src src-python">pd.crosstab(data.Smoker, data.Status, margins=<span style="color: #008b8b;">True</span>)
</pre>
</div>
<pre class="example">
Status Alive Dead All
Smoker
No 502 230 732
Yes 443 139 582
All 945 369 1314
</pre>
<p>
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!
</p>
<p>
Regardons maintenant les mêmes données sous forme de ratios par
rapport à <b>toute la population</b>. Par exemple, 10,6% des femmes de
l'étude étaient fumeuses avant de décéder.
</p>
<div class="org-src-container">
<pre class="src src-python">pd.crosstab(data.Smoker, data.Status, margins=<span style="color: #008b8b;">True</span>, normalize=<span style="color: #8b2252;">'all'</span>)*100
</pre>
</div>
<pre class="example">
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
</pre>
<p>
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%)
</p>
<p>
Déterminons maintenant les taux de mortalités pour chaque groupe
(fumeuse/ non fumeuse):
</p>
<div class="org-src-container">
<pre class="src src-python">pd.crosstab(data.Smoker, data.Status, margins=<span style="color: #008b8b;">True</span>, normalize=<span style="color: #8b2252;">'index'</span>)*100
</pre>
</div>
<pre class="example">
Status Alive Dead
Smoker
No 68.6 31.4
Yes 76.1 23.9
All 71.9 28.1
</pre>
<p>
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%)
</p>
<p>
Observons ce résultat surprenant sous forme graphique:
</p>
<div class="org-src-container">
<pre class="src src-python">pd.crosstab(data.Smoker, data.Status).plot.barh(stacked=<span style="color: #008b8b;">True</span>)
<span style="color: #a0522d;">figure_path</span> = <span style="color: #8b2252;">"data_as_stacked_bar_plot.png"</span>
plt.savefig(figure_path)
figure_path
</pre>
</div>
<div class="figure">
<p><img src="data_as_stacked_bar_plot.png" alt="data_as_stacked_bar_plot.png" />
</p>
</div>
<p>
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 <i>MELL-O-WELL</i> comme le recommande ce
médecin sur une publicité de 1930 ?!
</p>
<p>
<img src="https://www.vivelapub.fr/wp-content/uploads/2012/08/clop/176.jpg" alt="176.jpg" />
(<a href="https://tobacco.stanford.edu/tobacco_main/">source</a>)
</p>
</div>
</div>
<div id="outline-container-org3788f96" class="outline-3">
<h3 id="org3788f96"><span class="section-number-3">2.2</span> Partitionnement selon la classe d'âge</h3>
<div class="outline-text-3" id="text-2-2">
<p>
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.
</p>
<p>
Comme suggéré, groupons par tranches d'age les femmes de l'études:
</p>
<div class="org-src-container">
<pre class="src src-python"><span style="color: #b22222;"># </span><span style="color: #b22222;">On choisit 17 ans en borne inf&#233;rieur pour &#234;tre s&#251;r d'inclure le minimum (18)</span>
<span style="color: #a0522d;">data</span>[<span style="color: #8b2252;">'AgeCut'</span>] = pd.cut(data.Age, bins=[17, 34, 54, 64, data.Age.<span style="color: #483d8b;">max</span>()], labels=[<span style="color: #8b2252;">"18-34 ans"</span>, <span style="color: #8b2252;">"34-54 ans"</span>, <span style="color: #8b2252;">"55-64 ans"</span>, <span style="color: #8b2252;">"+ de 65 ans"</span>])
data.AgeCut
</pre>
</div>
<pre class="example">
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 &lt; 34-54 ans &lt; 55-64 ans &lt; + de 65 ans]
</pre>
<p>
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:
</p>
<div class="org-src-container">
<pre class="src src-python"><span style="color: #a020f0;">print</span>(<span style="color: #8b2252;">"Avec Quartiles:"</span>, *pd.qcut(data.Age, 4, retbins=<span style="color: #008b8b;">True</span>)[1].astype(<span style="color: #483d8b;">int</span>), sep=<span style="color: #8b2252;">'|'</span>)
<span style="color: #a020f0;">print</span>(<span style="color: #8b2252;">"Sans Quartiles:"</span>, *[18, 34, 54, 64, 89], sep=<span style="color: #8b2252;">'|'</span>)
</pre>
</div>
<pre class="example">
Avec Quartiles:|18|31|44|60|89
Sans Quartiles:|18|34|54|64|89
</pre>
<p>
Dénombrons les femmes ayant survécu/décédé selon leur tranche d'age
et si elles fumaient ou non:
</p>
<div class="org-src-container">
<pre class="src src-python">pd.crosstab([data.Smoker, data.AgeCut], data.Status, margins=<span style="color: #008b8b;">True</span>)
</pre>
</div>
<pre class="example">
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
</pre>
<p>
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%) <b>lors de la première étude</b>.
</p>
<p>
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é ?!
</p>
<div class="org-src-container">
<pre class="src src-python">pd.crosstab([data.AgeCut, data.Smoker], data.Status, margins=<span style="color: #008b8b;">True</span>, normalize=<span style="color: #8b2252;">'all'</span>)*100
</pre>
</div>
<pre class="example">
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
</pre>
<p>
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
<b>aussi pour les plus agés</b>!(~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 ?
</p>
<div class="org-src-container">
<pre class="src src-python">pd.crosstab([data.AgeCut, data.Smoker], data.Status, margins=<span style="color: #008b8b;">True</span>, normalize=<span style="color: #8b2252;">'index'</span>)*100
</pre>
</div>
<pre class="example">
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
</pre>
<p>
Visuallement, c'est plus parlant:
</p>
<div class="org-src-container">
<pre class="src src-python">pd.crosstab([data.AgeCut, data.Smoker], data.Status).plot.barh(stacked=<span style="color: #008b8b;">True</span>)
<span style="color: #a0522d;">figure_path</span> = <span style="color: #8b2252;">"data_as_stacked_bar_plot_agecut.png"</span>
plt.savefig(figure_path, bbox_inches=<span style="color: #8b2252;">'tight'</span>)
figure_path
</pre>
</div>
<p>
<img src="data_as_stacked_bar_plot_agecut.png" alt="data_as_stacked_bar_plot_agecut.png" />
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.
</p>
<p>
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.
</p>
<p>
Au Royaume-Uni, l'espérance de vie à la naissance en 1974 était de 72,5 ans
(<a href="https://donnees.banquemondiale.org/indicateur/SP.DYN.LE00.IN?locations=GB">source</a>). 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.
</p>
</div>
</div>
<div id="outline-container-org25b0d41" class="outline-3">
<h3 id="org25b0d41"><span class="section-number-3">2.3</span> Modélisation de la probabilité de décès en fonction de l'âge</h3>
<div class="outline-text-3" id="text-2-3">
<p>
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).
</p>
<p>
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).
</p>
<p>
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 ?
</p>
<div class="org-src-container">
<pre class="src src-python"><span style="color: #a020f0;">import</span> seaborn <span style="color: #a020f0;">as</span> sns
plt.figure(figsize=(9, 7))
sns.violinplot(x=<span style="color: #8b2252;">'Smoker'</span>, y=<span style="color: #8b2252;">'Age'</span>, hue=<span style="color: #8b2252;">'Status'</span>, data=data, split=<span style="color: #008b8b;">True</span>, inner=<span style="color: #8b2252;">'quartile'</span>)
sns.despine()
<span style="color: #a0522d;">figure_path</span> = <span style="color: #8b2252;">"smoker_survival_violin.png"</span>
plt.savefig(figure_path, bbox_inches=<span style="color: #8b2252;">'tight'</span>)
figure_path
</pre>
</div>
<div class="figure">
<p><img src="smoker_survival_violin.png" alt="smoker_survival_violin.png" />
</p>
</div>
<p>
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:
</p>
<div class="org-src-container">
<pre class="src src-python"><span style="color: #b22222;"># </span><span style="color: #b22222;">On d&#233;finit une variable Death(0 si survie, 1 si d&#233;c&#232;s)</span>
<span style="color: #a0522d;">data</span>[<span style="color: #8b2252;">'Death'</span>] = data.Status.cat.codes
data[[<span style="color: #8b2252;">'Status'</span>, <span style="color: #8b2252;">'Death'</span>]]
</pre>
</div>
<pre class="example">
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]
</pre>
<p>
Observons rapidement ce qu'une modélisation avec une régression
logistique sur chaque groupe fumeuse/non fumeuse pourrait donner:
</p>
<div class="org-src-container">
<pre class="src src-python">plt.close()
plt.figure(figsize=(12, 15))
sns.lmplot(x=<span style="color: #8b2252;">"Age"</span>, y=<span style="color: #8b2252;">"Death"</span>, data=data, logistic=<span style="color: #008b8b;">True</span>, hue=<span style="color: #8b2252;">'Smoker'</span>, y_jitter=.03, aspect=1.5, markers=[<span style="color: #8b2252;">'o'</span>, <span style="color: #8b2252;">'x'</span>]);
sns.despine()
<span style="color: #a0522d;">figure_path</span> = <span style="color: #8b2252;">"smoker_survival_logistic.png"</span>
plt.savefig(figure_path, bbox_inches=<span style="color: #8b2252;">'tight'</span>)
figure_path
</pre>
</div>
<div class="figure">
<p><img src="smoker_survival_logistic.png" alt="smoker_survival_logistic.png" />
</p>
</div>
<p>
Plusieurs remarques intéressantes:
</p>
<ul class="org-ul">
<li>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.</li>
<li>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</li>
<li>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.</li>
</ul>
<p>
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
<a href="https://github.com/mwaskom/seaborn/issues/655#issuecomment-125073496">ne
sera jamais permis!</a>
</p>
<p>
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:
</p>
<div class="org-src-container">
<pre class="src src-python"><span style="color: #a020f0;">from</span> sklearn.model_selection <span style="color: #a020f0;">import</span> cross_val_score, cross_validate
<span style="color: #a020f0;">from</span> sklearn.model_selection <span style="color: #a020f0;">import</span> RepeatedStratifiedKFold
<span style="color: #a020f0;">def</span> <span style="color: #0000ff;">train</span>(model, model_name, data):
<span style="color: #8b2252;">"""Entraine un mod&#232;le pour chaque partition(Fumeuse/Non Fumeuse).</span>
<span style="color: #8b2252;"> Ajoute &#233;galement 2 colonnes &#224; data contenant les pr&#233;dictions du mod&#232;les ainsi que les probabilit&#233;s"""</span>
<span style="color: #b22222;"># </span><span style="color: #b22222;">On stocke les classifieurs pour chaque groupe (fumeuse/non fumeuse)</span>
<span style="color: #a0522d;">clf</span> = <span style="color: #483d8b;">dict</span>(Yes=<span style="color: #008b8b;">None</span>, No=<span style="color: #008b8b;">None</span>, name=model_name)
<span style="color: #a020f0;">for</span> is_smoker <span style="color: #a020f0;">in</span> [<span style="color: #8b2252;">'Yes'</span>, <span style="color: #8b2252;">'No'</span>]:
<span style="color: #b22222;"># </span><span style="color: #b22222;">D&#233;finition des jeux de donn&#233;es pour la validation crois&#233;e</span>
<span style="color: #a0522d;">cv</span> = RepeatedStratifiedKFold(n_splits=6, n_repeats=3, random_state=1)
<span style="color: #a0522d;">X</span> = data[data.Smoker == is_smoker][<span style="color: #8b2252;">'Age'</span>].to_frame()
<span style="color: #a0522d;">targets</span> = data[data.Smoker == is_smoker][<span style="color: #8b2252;">'Death'</span>]
<span style="color: #b22222;"># </span><span style="color: #b22222;">Entrainement</span>
<span style="color: #a0522d;">scores</span> = cross_validate(model, X, targets, scoring=<span style="color: #8b2252;">'roc_auc'</span>, cv=cv, n_jobs=-1, return_estimator=<span style="color: #008b8b;">True</span>)
<span style="color: #b22222;"># </span><span style="color: #b22222;">On r&#233;cup&#232;re le meilleur classifieur</span>
clf.update({is_smoker: scores[<span style="color: #8b2252;">'estimator'</span>][scores[<span style="color: #8b2252;">'test_score'</span>].argmax()]})
<span style="color: #b22222;"># </span><span style="color: #b22222;">Ajout des pr&#233;dictions et probabilit&#233;s du mod&#232;le entrain&#233; pour la partition en cours(fumeuse ou non fumeuse)</span>
<span style="color: #b22222;"># </span><span style="color: #b22222;">Les probabilit&#233;s retourn&#233;es sont [P(Death = 0), P(Death = 1)]. On stocke donc le deuxi&#232;me &#233;l&#233;ment.</span>
<span style="color: #b22222;"># </span><span style="color: #b22222;">Visualiser ordre des classes{0,1} dans clf[is_smoker].classes_ pour retrouver cette information</span>
<span style="color: #a0522d;">data.loc</span>[data.Smoker == is_smoker, f<span style="color: #8b2252;">'DeathProba{model_name}'</span>] = clf[is_smoker].predict_proba(data[data.Smoker == is_smoker][<span style="color: #8b2252;">'Age'</span>].to_frame())[:,1]
<span style="color: #a020f0;">return</span> clf
</pre>
</div>
<p>
Nous allons également entraîner un modèle d'arbre de décision, qui
peut se retrouver assez performant pour la classification binaire.
</p>
<div class="org-src-container">
<pre class="src src-python"><span style="color: #a020f0;">from</span> sklearn.linear_model <span style="color: #a020f0;">import</span> LogisticRegression
<span style="color: #a020f0;">from</span> sklearn.tree <span style="color: #a020f0;">import</span> DecisionTreeClassifier
<span style="color: #b22222;"># </span><span style="color: #b22222;">train models</span>
<span style="color: #a0522d;">model</span> = LogisticRegression(solver=<span style="color: #8b2252;">'lbfgs'</span>, class_weight=<span style="color: #8b2252;">'balanced'</span>, C=0.1)
<span style="color: #a0522d;">clf_LR</span> = train(model, <span style="color: #8b2252;">"LR"</span>, data)
<span style="color: #a0522d;">model</span> = DecisionTreeClassifier(class_weight=<span style="color: #8b2252;">'balanced'</span>, min_samples_leaf=5)
<span style="color: #a0522d;">clf_tree</span> = train(model, <span style="color: #8b2252;">"Tree"</span>, data)
</pre>
</div>
<p>
Observons les résultats
</p>
<div class="org-src-container">
<pre class="src src-python"><span style="color: #a020f0;">import</span> numpy <span style="color: #a020f0;">as</span> np
<span style="color: #a020f0;">from</span> scipy.special <span style="color: #a020f0;">import</span> expit
plt.clf()
<span style="color: #a0522d;">colors</span> = <span style="color: #483d8b;">dict</span>(Yes=<span style="color: #8b2252;">'green'</span>, No=<span style="color: #8b2252;">'Red'</span>)
<span style="color: #a020f0;">for</span> is_smoker <span style="color: #a020f0;">in</span> [<span style="color: #8b2252;">'Yes'</span>, <span style="color: #8b2252;">'No'</span>]:
<span style="color: #a0522d;">clf</span> = clf_LR[is_smoker]
<span style="color: #a0522d;">age</span>, <span style="color: #a0522d;">death</span> = data[data.Smoker == is_smoker][[<span style="color: #8b2252;">'Age'</span>, <span style="color: #8b2252;">'Death'</span>]].T.values
<span style="color: #b22222;"># </span><span style="color: #b22222;">Prepare for drawing logistic curve estimate</span>
<span style="color: #a0522d;">age_values</span> = np.linspace(-0, 100, 300)
<span style="color: #a0522d;">logit</span> = expit(age_values * clf.coef_ + clf.intercept_).ravel()
<span style="color: #b22222;"># </span><span style="color: #b22222;">Draw</span>
plt.scatter(age, death, color=colors[is_smoker], marker=<span style="color: #8b2252;">'+'</span>)
plt.plot(age_values, logit, color=colors[is_smoker], linewidth=2)
<span style="color: #a0522d;">figure_path</span> = <span style="color: #8b2252;">'smoker_survival_logistic_custom.png'</span>
plt.savefig(figure_path)
figure_path
</pre>
</div>
<p>
<img src="smoker_survival_logistic_custom.png" alt="smoker_survival_logistic_custom.png" />
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).
</p>
<p>
Regardons quelques métriques usuelles
</p>
<div class="org-src-container">
<pre class="src src-python"><span style="color: #a020f0;">from</span> sklearn <span style="color: #a020f0;">import</span> metrics
<span style="color: #a0522d;">metrics_dict</span> = <span style="color: #483d8b;">dict</span>(accuracy=metrics.accuracy_score,
precision=metrics.precision_score,
recall=metrics.recall_score)
<span style="color: #a020f0;">for</span> clf <span style="color: #a020f0;">in</span> [clf_LR, clf_tree]:
<span style="color: #a020f0;">print</span>(f<span style="color: #8b2252;">"&lt;&lt;&lt;clf is {clf['Yes'].__class__.__name__}&gt;&gt;&gt;"</span>)
<span style="color: #a020f0;">for</span> is_smoker <span style="color: #a020f0;">in</span> [<span style="color: #8b2252;">'Yes'</span>, <span style="color: #8b2252;">'No'</span>]:
<span style="color: #a020f0;">print</span>(f<span style="color: #8b2252;">"------Smoker is {is_smoker}------"</span>)
<span style="color: #a0522d;">y_true</span> = data[data.Smoker == is_smoker][<span style="color: #8b2252;">'Death'</span>]
<span style="color: #a0522d;">y_pred</span> = clf[is_smoker].predict(data[data.Smoker == is_smoker][<span style="color: #8b2252;">'Age'</span>].to_frame())
<span style="color: #a020f0;">for</span> metric,score <span style="color: #a020f0;">in</span> metrics_dict.items():
<span style="color: #a020f0;">print</span>(f<span style="color: #8b2252;">"{metric:&gt;10} : {score(y_true, y_pred):.3}"</span>)
<span style="color: #b22222;">#</span><span style="color: #b22222;">print(metrics.confusion_matrix(y_true, y_pred)))</span>
</pre>
</div>
<pre class="example">
&lt;&lt;&lt;clf is LogisticRegression&gt;&gt;&gt;
------Smoker is Yes------
accuracy : 0.722
precision : 0.448
recall : 0.712
------Smoker is No------
accuracy : 0.832
precision : 0.683
recall : 0.87
&lt;&lt;&lt;clf is DecisionTreeClassifier&gt;&gt;&gt;
------Smoker is Yes------
accuracy : 0.789
precision : 0.537
recall : 0.827
------Smoker is No------
accuracy : 0.851
precision : 0.711
recall : 0.887
</pre>
<div class="org-src-container">
<pre class="src src-python"><span style="color: #a020f0;">from</span> itertools <span style="color: #a020f0;">import</span> product
plt.clf()
<span style="color: #a020f0;">for</span> clf,is_smoker <span style="color: #a020f0;">in</span> product([clf_LR, clf_tree], [<span style="color: #8b2252;">'Yes'</span>, <span style="color: #8b2252;">'No'</span>]):
<span style="color: #a0522d;">y_true</span> = data[data.Smoker == is_smoker][<span style="color: #8b2252;">'Death'</span>]
<span style="color: #a0522d;">y_pred_proba</span> = data[data.Smoker == is_smoker][f<span style="color: #8b2252;">"DeathProba{clf['name']}"</span>]
<span style="color: #a0522d;">fpr</span>, <span style="color: #a0522d;">tpr</span>, <span style="color: #a0522d;">_</span> = metrics.roc_curve(y_true, y_pred_proba)
<span style="color: #a0522d;">auc</span> = metrics.roc_auc_score(y_true, y_pred_proba)
plt.plot(fpr,tpr,label=f<span style="color: #8b2252;">"{clf['name']}(Smoker={is_smoker}), auc={auc:.2}"</span>)
plt.legend(loc=4)
plt.xlabel(<span style="color: #8b2252;">'Taux de faux positifs'</span>)
plt.ylabel(<span style="color: #8b2252;">'Taux de vrais positifs'</span>)
plt.title(<span style="color: #8b2252;">'Courbe ROC'</span>)
<span style="color: #a0522d;">figure_path</span> = <span style="color: #8b2252;">"roc.png"</span>
plt.savefig(figure_path, bbox_inches=<span style="color: #8b2252;">'tight'</span>);
plt.close()
<span style="color: #a020f0;">return</span> figure_path
</pre>
</div>
<div class="figure">
<p><img src="roc.png" alt="roc.png" />
</p>
</div>
<p>
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.
</p>
<p>
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.
</p>
</div>
</div>
</div>
</div>
<div id="postamble" class="status">
<p class="date">Date: 2020-04-21</p>
<p class="author">Auteur: Guillaume</p>
<p class="date">Created: 2020-04-24 ven. 09:47</p>
<p class="validation"><a href="http://validator.w3.org/check?uri=referer">Validate</a></p>
</div>
</body>
</html>
#+TITLE: Sujet 6: Autour du Paradoxe de Simpson
#+AUTHOR: Guillaume
#+DATE: 2020-04-21
#+LANGUAGE: fr
# #+PROPERTY: header-args :session **python** :exports both
#+PROPERTY: header-args :session *python* :exports both
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="http://www.pirilampo.org/styles/readtheorg/css/htmlize.css"/>
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="http://www.pirilampo.org/styles/readtheorg/css/readtheorg.css"/>
#+HTML_HEAD: <script src="https://ajax.googleapis.com/ajax/libs/jquery/2.1.3/jquery.min.js"></script>
#+HTML_HEAD: <script src="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.4/js/bootstrap.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="http://www.pirilampo.org/styles/lib/js/jquery.stickytableheaders.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="http://www.pirilampo.org/styles/readtheorg/js/readtheorg.js"></script>
#+OPTIONS: H:4 toc:3 p:t
* 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.
([[https://gitlab.inria.fr/learninglab/mooc-rr/mooc-rr-ressources/blob/master//module3/Practical_session/Subject6_fr.md][source]])
* Mission
** 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
#+begin_src python :results value :exports both
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
#+end_src
#+RESULTS:
#+begin_example
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]
#+end_example
#+begin_src python :results value :exports both
data.describe(include='all')
#+end_src
#+RESULTS:
#+begin_example
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
#+end_example
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:
#+begin_src python :results value :exports both
pd.crosstab(data.Smoker, data.Status, margins=True)
#+end_src
#+RESULTS:
: 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.
#+begin_src python :results value :exports both
pd.crosstab(data.Smoker, data.Status, margins=True, normalize='all')*100
#+end_src
#+RESULTS:
: 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):
#+begin_src python :results value :exports both
pd.crosstab(data.Smoker, data.Status, margins=True, normalize='index')*100
#+end_src
#+RESULTS:
: 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:
#+begin_src python :results file :exports both
pd.crosstab(data.Smoker, data.Status).plot.barh(stacked=True)
figure_path = "data_as_stacked_bar_plot.png"
plt.savefig(figure_path)
figure_path
#+end_src
#+RESULTS:
[[file:data_as_stacked_bar_plot.png]]
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 ?!
[[https://www.vivelapub.fr/wp-content/uploads/2012/08/clop/176.jpg]]
([[https://tobacco.stanford.edu/tobacco_main/][source]])
** 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:
#+begin_src python :results value :exports both
# 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
#+end_src
#+RESULTS:
#+begin_example
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]
#+end_example
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:
#+begin_src python :results output :exports both
print("Avec Quartiles:", *pd.qcut(data.Age, 4, retbins=True)[1].astype(int), sep='|')
print("Sans Quartiles:", *[18, 34, 54, 64, 89], sep='|')
#+end_src
#+RESULTS:
: 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:
#+begin_src python :results value :exports both
pd.crosstab([data.Smoker, data.AgeCut], data.Status, margins=True)
#+end_src
#+RESULTS:
#+begin_example
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
#+end_example
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é ?!
#+begin_src python :results value :exports both
pd.crosstab([data.AgeCut, data.Smoker], data.Status, margins=True, normalize='all')*100
#+end_src
#+RESULTS:
#+begin_example
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
#+end_example
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 ?
#+begin_src python :results value :exports both
pd.crosstab([data.AgeCut, data.Smoker], data.Status, margins=True, normalize='index')*100
#+end_src
#+RESULTS:
#+begin_example
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
#+end_example
Visuallement, c'est plus parlant:
#+begin_src python :results file :exports both
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
#+end_src
#+RESULTS:
[[file:data_as_stacked_bar_plot_agecut.png]]
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
([[https://donnees.banquemondiale.org/indicateur/SP.DYN.LE00.IN?locations=GB][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.
** 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 ?
#+begin_src python :results file :exports both
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
#+end_src
#+RESULTS:
[[file:smoker_survival_violin.png]]
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:
#+begin_src python :results value :exports both
# On définit une variable Death(0 si survie, 1 si décès)
data['Death'] = data.Status.cat.codes
data[['Status', 'Death']]
#+end_src
#+RESULTS:
#+begin_example
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]
#+end_example
Observons rapidement ce qu'une modélisation avec une régression
logistique sur chaque groupe fumeuse/non fumeuse pourrait donner:
#+begin_src python :results file :exports both
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
#+end_src
#+RESULTS:
[[file:smoker_survival_logistic.png]]
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
[[https://github.com/mwaskom/seaborn/issues/655#issuecomment-125073496][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:
#+begin_src python :results silent :exports both
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
#+end_src
#+RESULTS:
Nous allons également entraîner un modèle d'arbre de décision, qui
peut se retrouver assez performant pour la classification binaire.
#+begin_src python :results output :exports both
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)
#+end_src
#+RESULTS:
Observons les résultats
#+begin_src python :results silent :exports none
"""from sklearn.tree import plot_tree
plt.close()
plt.figure(figsize=(4,3), dpi=400)
plot_tree(clf_tree['No'], filled=True, feature_names=['Age'])
plt.tight_layout()
plt.show()"""
#+end_src
#+RESULTS:
#+begin_src python :results file :exports both
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
#+end_src
#+RESULTS:
[[file:smoker_survival_logistic_custom.png]]
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
#+RESULTS:
#+begin_src python :results output :exports both
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)))
#+end_src
#+RESULTS:
#+begin_example
<<<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
#+end_example
#+begin_src python :results silent :exports code
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
#+end_src
#+RESULTS:
[[file:roc.png]]
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.
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment