Commit 85b076e2 authored by Jade Bolaty's avatar Jade Bolaty

sujet 6 paradoxe simpson orgmode r avec tout ce qui va avec

parent d460d83e
#+TITLE: À propos du calcul de $\pi$
#+AUTHOR: Arnaud Legrand
#+DATE: <2024-11-02 sam.>
#+LANGUAGE: fr
# #+PROPERTY: header-args :eval never-export
#+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: num:nil
* Table des matières
+ [[uno][1 En demandant à la lib maths]]
+ [[dos][2 En utilisant la méthode des aiguilles de Buffon]]
+ [[tres][3 Avec un argument "fréquentiel" de surface]]
* <<uno>>1 En demandant à la lib maths
Mon ordinateur m'indique que $\pi$ vaut /approximativement/
#+begin_src R :results output :session *R* :exports both
pi
#+end_src
#+RESULTS:
: [1] 3.141593
* <<dos>>2 En utilisant la méthode des aiguilles de Buffon
Mais calculé avec la *méthode* des [[https://fr.wikipedia.org/wiki/Aiguille_de_Buffon][aiguilles de Buffon]], on obtiendrait
comme *approximation* :
#+begin_src R :results output :session *R* :exports both
set.seed(42)
N = 100000
x = runif(N)
theta = pi/2*runif(N)
2/(mean(x+sin(theta)>1))
#+end_src
#+RESULTS:
: [1] 3.14327
* <<tres>>3 Avec un argument "fréquentiel" de surface
Sinon, une méthode plus simple à comprendre et ne faisant pas
intervenir d'appel à lz fonction sinus se base sur le fait que si $X \sim
U(0,1)$ et $Y \sim U(0,1)$ alors $P[X² + Y² <= 1] = \pi/4$ (voir [[https://fr.wikipedia.org/wiki/M%C3%A9thode_de_Monte-Carlo#D%C3%A9termination_de_la_valeur_de_%CF%80][méthode de
Monte Carlo sur Wikipedia]]). Le code suivant illustre ce fait :
#+begin_src R :results file graphics :file (org-babel-temp-file "figure" ".png") :exports both :width 600 :height 400 :session *R*
set.seed(42)
N = 1000
df = data.frame(X = runif(N), Y = runif(N))
df$Accept = (df$X**2 + df$Y**2 <=1)
library(ggplot2)
ggplot(df, aes(x=X,y=Y,color=Accept)) + geom_point(alpha=.2) + coord_fixed() + theme_bw()
#+end_src
#+RESULTS:
[[file:c:/Users/Jade/AppData/Local/Temp/babel-v0Y8bg/figurez1fDhU.png]]
Il est alors aisé d'obtenir une approximation (pas terrible) de $\pi$ en
comptant combien de fois, en moyenne, $X^2 + Y^2$ est inférieur à 1 :
#+begin_src R :results output :session *R* :exports both
4*mean(df$Accept)
#+end_src
#+RESULTS:
: [1] 3.156
Jade@LORDI-RKHGKLV94.17136:1730227632
\ No newline at end of file
Jade@LORDI-RKHGKLV94.15420:1730565465
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
#+TITLE: Votre titre
#+AUTHOR: Votre nom
#+DATE: La date du jour
#+LANGUAGE: fr
# #+PROPERTY: header-args :eval never-export
#+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>
#+begin_src R :results output :session *R* :exports both
mylist <- c(14.0, 7.6, 11.2, 12.8, 12.5, 9.9, 14.9, 9.4, 16.9, 10.2, 14.9, 18.1, 7.3, 9.8, 10.9,12.2, 9.9, 2.9, 2.8, 15.4, 15.7, 9.7, 13.1, 13.2, 12.3, 11.7, 16.0, 12.4, 17.9, 12.2, 16.2, 18.7, 8.9, 11.9, 12.1, 14.6, 12.1, 4.7, 3.9, 16.9, 16.8, 11.3, 14.4, 15.7, 14.0, 13.6, 18.0, 13.6, 19.9, 13.7, 17.0, 20.5, 9.9, 12.5, 13.2, 16.1, 13.5, 6.3, 6.4, 17.6, 19.1, 12.8, 15.5, 16.3, 15.2, 14.6, 19.1, 14.4, 21.4, 15.1, 19.6, 21.7, 11.3, 15.0, 14.3, 16.8, 14.0, 6.8, 8.2, 19.9, 20.4, 14.6, 16.4, 18.7, 16.8, 15.8, 20.4, 15.8, 22.4, 16.2, 20.3, 23.4, 12.1, 15.5, 15.4, 18.4, 15.7, 10.2, 8.9, 21.0)
mean(mylist)
min(mylist)
max(mylist)
median(mylist)
sd(mylist)
#+end_src
#+RESULTS:
: [1] 14.113
: [1] 2.8
: [1] 23.4
: [1] 14.5
: [1] 4.334094
#+TITLE: Votre titre
#+AUTHOR: Votre nom
#+DATE: La date du jour
#+LANGUAGE: fr
# #+PROPERTY: header-args :eval never-export
#+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>
#+begin_src python :results output :session :exports both
import numpy as np
mylist = [14.0, 7.6, 11.2, 12.8, 12.5, 9.9, 14.9, 9.4, 16.9, 10.2, 14.9, 18.1, 7.3, 9.8, 10.9,12.2, 9.9, 2.9, 2.8, 15.4, 15.7, 9.7, 13.1, 13.2, 12.3, 11.7, 16.0, 12.4, 17.9, 12.2, 16.2, 18.7, 8.9, 11.9, 12.1, 14.6, 12.1, 4.7, 3.9, 16.9, 16.8, 11.3, 14.4, 15.7, 14.0, 13.6, 18.0, 13.6, 19.9, 13.7, 17.0, 20.5, 9.9, 12.5, 13.2, 16.1, 13.5, 6.3, 6.4, 17.6, 19.1, 12.8, 15.5, 16.3, 15.2, 14.6, 19.1, 14.4, 21.4, 15.1, 19.6, 21.7, 11.3, 15.0, 14.3, 16.8, 14.0, 6.8, 8.2, 19.9, 20.4, 14.6, 16.4, 18.7, 16.8, 15.8, 20.4, 15.8, 22.4, 16.2, 20.3, 23.4, 12.1, 15.5, 15.4, 18.4, 15.7, 10.2, 8.9, 21.0]
print(np.mean(mylist))
print(np.min(mylist))
print(np.max(mylist))
print(np.median(mylist))
print(np.std(mylist,ddof=1))
#+end_src
#+RESULTS:
: 14.113000000000001
: 2.8
: 23.4
: 14.5
: 4.334094455301447
Jade@LORDI-RKHGKLV94.16232:1730565748
\ No newline at end of file
Jade@LORDI-RKHGKLV94.16232:1730565748
\ No newline at end of file
#+TITLE: Analyse du risque de défaillance des joints toriques de la navette Challenger
#+AUTHOR: Arnaud Legrand
#+LANGUAGE: fr
#+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>
#+LATEX_HEADER: \usepackage[utf8]{inputenc}
#+LATEX_HEADER: \usepackage[T1]{fontenc}
#+LATEX_HEADER: \usepackage[a4paper,margin=.8in]{geometry}
#+LATEX_HEADER: \usepackage[french]{babel}
# #+PROPERTY: header-args :session :exports both
Le 27 Janvier 1986, veille du décollage de la navette /Challenger/, eu
lieu une télé-conférence de trois heures entre les ingénieurs de la
Morton Thiokol (constructeur d'un des moteurs) et de la NASA. La
discussion portait principalement sur les conséquences de la
température prévue au moment du décollage de 31°F (juste en dessous de
0°C) sur le succès du vol et en particulier sur la performance des
joints toriques utilisés dans les moteurs. En effet, aucun test
n'avait été effectué à cette température.
L'étude qui suit reprend donc une partie des analyses effectuées cette
nuit là et dont l'objectif était d'évaluer l'influence potentielle de
la température et de la pression à laquelle sont soumis les joints
toriques sur leur probabilité de dysfonctionnement. Pour cela, nous
disposons des résultats des expériences réalisées par les ingénieurs
de la NASA durant les 6 années précédant le lancement de la navette
Challenger.
* Chargement des données
Nous commençons donc par charger ces données:
#+begin_src R :results output :session *R* :exports both
data = read.csv("C:/Users/Jade/mooc-rr/module2/exo5/shuttle.csv",header=T)
data
#+end_src
#+RESULTS:
#+begin_example
Date Count Temperature Pressure Malfunction
1 4/12/81 6 66 50 0
2 11/12/81 6 70 50 1
3 3/22/82 6 69 50 0
4 11/11/82 6 68 50 0
5 4/04/83 6 67 50 0
6 6/18/82 6 72 50 0
7 8/30/83 6 73 100 0
8 11/28/83 6 70 100 0
9 2/03/84 6 57 200 1
10 4/06/84 6 63 200 1
11 8/30/84 6 70 200 1
12 10/05/84 6 78 200 0
13 11/08/84 6 67 200 0
14 1/24/85 6 53 200 2
15 4/12/85 6 67 200 0
16 4/29/85 6 75 200 0
17 6/17/85 6 70 200 0
18 7/29/85 6 81 200 0
19 8/27/85 6 76 200 0
20 10/03/85 6 79 200 0
21 10/30/85 6 75 200 2
22 11/26/85 6 76 200 0
23 1/12/86 6 58 200 1
#+end_example
Le jeu de données nous indique la date de l'essai, le nombre de joints
toriques mesurés (il y en a 6 sur le lançeur principal), la
température (en Fahrenheit) et la pression (en psi), et enfin le
nombre de dysfonctionnements relevés.
* Inspection graphique des données
Les vols où aucun incident n'est relevé n'apportant aucune information
sur l'influence de la température ou de la pression sur les
dysfonctionnements, nous nous concentrons sur les expériences où au
moins un joint a été défectueux.
#+begin_src R :results output :session *R* :exports both
data = data[data$Malfunction>0,]
data
#+end_src
#+RESULTS:
: Date Count Temperature Pressure Malfunction
: 2 11/12/81 6 70 50 1
: 9 2/03/84 6 57 200 1
: 10 4/06/84 6 63 200 1
: 11 8/30/84 6 70 200 1
: 14 1/24/85 6 53 200 2
: 21 10/30/85 6 75 200 2
: 23 1/12/86 6 58 200 1
Très bien, nous avons une variabilité de température importante mais
la pression est quasiment toujours égale à 200, ce qui devrait
simplifier l'analyse.
Comment la fréquence d'échecs varie-t-elle avec la température ?
#+begin_src R :results file graphics :file "freq_temp.png" :exports both :width 600 :height 400 :session *R*
plot(data=data, Malfunction/Count ~ Temperature, ylim=c(0,1))
#+end_src
#+RESULTS:
[[file:freq_temp.png]]
À première vue, ce n'est pas flagrant mais bon, essayons quand même
d'estimer l'impact de la température $t$ sur la probabilité de
dysfonctionnements d'un joint.
* Estimation de l'influence de la température
Supposons que chacun des 6 joints toriques est endommagé avec la même
probabilité et indépendamment des autres et que cette probabilité ne
dépend que de la température. Si on note $p(t)$ cette probabilité, le
nombre de joints $D$ dysfonctionnant lorsque l'on effectue le vol à
température $t$ suit une loi binomiale de paramètre $n=6$ et
$p=p(t)$. Pour relier $p(t)$ à $t$, on va donc effectuer une
régression logistique.
#+begin_src R :results output :session *R* :exports both
logistic_reg = glm(data=data, Malfunction/Count ~ Temperature, weights=Count,
family=binomial(link='logit'))
summary(logistic_reg)
#+end_src
#+RESULTS:
#+begin_example
Call:
glm(formula = Malfunction/Count ~ Temperature, family = binomial(link = "logit"),
data = data, weights = Count)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.389528 3.195752 -0.435 0.664
Temperature 0.001416 0.049773 0.028 0.977
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1.3347 on 6 degrees of freedom
Residual deviance: 1.3339 on 5 degrees of freedom
AIC: 18.894
Number of Fisher Scoring iterations: 4
#+end_example
L'estimateur le plus probable du paramètre de température est 0.001416
et l'erreur standard de cet estimateur est de 0.049, autrement dit on
ne peut pas distinguer d'impact particulier et il faut prendre nos
estimations avec des pincettes.
* Estimation de la probabilité de dysfonctionnant des joints toriques
La température prévue le jour du décollage est de 31°F. Essayons
d'estimer la probabilité de dysfonctionnement des joints toriques à
cette température à partir du modèle que nous venons de construire:
#+begin_src R :results file graphics :file "proba_estimate.png" :exports both :width 600 :height 400 :session *R*
# shuttle=shuttle[shuttle$r!=0,]
tempv = seq(from=30, to=90, by = .5)
rmv <- predict(logistic_reg,list(Temperature=tempv),type="response")
plot(tempv,rmv,type="l",ylim=c(0,1))
points(data=data, Malfunction/Count ~ Temperature)
#+end_src
#+RESULTS:
[[file:proba_estimate.png]]
Comme on pouvait s'attendre au vu des données initiales, la
température n'a pas d'impact notable sur la probabilité d'échec des
joints toriques. Elle sera d'environ 0.2, comme dans les essais
précédents où nous il y a eu défaillance d'au moins un joint. Revenons
à l'ensemble des données initiales pour estimer la probabilité de
défaillance d'un joint:
#+begin_src R :results output :session *R* :exports both
data_full = read.csv("C:/Users/Jade/mooc-rr/module2/exo5/shuttle.csv",header=T)
sum(data_full$Malfunction)/sum(data_full$Count)
#+end_src
#+RESULTS:
: [1] 0.06521739
Cette probabilité est donc d'environ $p=0.065$, sachant qu'il existe
un joint primaire un joint secondaire sur chacune des trois parties du
lançeur, la probabilité de défaillance des deux joints d'un lançeur
est de $p^2 \approx 0.00425$. La probabilité de défaillance d'un des
lançeur est donc de $1-(1-p^2)^3 \approx 1.2%$. Ça serait vraiment
pas de chance... Tout est sous contrôle, le décollage peut donc avoir
lieu demain comme prévu.
Seulement, le lendemain, la navette Challenger explosera et emportera
avec elle ses sept membres d'équipages. L'opinion publique est
fortement touchée et lors de l'enquête qui suivra, la fiabilité des
joints toriques sera directement mise en cause. Au delà des problèmes
de communication interne à la NASA qui sont pour beaucoup dans ce
fiasco, l'analyse précédente comporte (au moins) un petit
problème... Saurez-vous le trouver ? Vous êtes libre de modifier cette
analyse et de regarder ce jeu de données sous tous les angles afin
d'expliquer ce qui ne va pas.
Jade@LORDI-RKHGKLV94.11112:1731010242
\ No newline at end of file
This diff is collapsed.
......@@ -36,36 +36,36 @@ Challenger.
* Chargement des données
Nous commençons donc par charger ces données:
#+begin_src R :results output :session *R* :exports both
data = read.csv("shuttle.csv",header=T)
data = read.csv("C:/Users/Jade/mooc-rr/module2/exo5/shuttle.csv",header=T)
data
#+end_src
#+RESULTS:
#+begin_example
Date Count Temperature Pressure Malfunction
1 4/12/81 6 66 50 0
2 11/12/81 6 70 50 1
3 3/22/82 6 69 50 0
4 11/11/82 6 68 50 0
5 4/04/83 6 67 50 0
6 6/18/82 6 72 50 0
7 8/30/83 6 73 100 0
8 11/28/83 6 70 100 0
9 2/03/84 6 57 200 1
10 4/06/84 6 63 200 1
11 8/30/84 6 70 200 1
12 10/05/84 6 78 200 0
13 11/08/84 6 67 200 0
14 1/24/85 6 53 200 2
15 4/12/85 6 67 200 0
16 4/29/85 6 75 200 0
17 6/17/85 6 70 200 0
18 7/29/85 6 81 200 0
19 8/27/85 6 76 200 0
20 10/03/85 6 79 200 0
21 10/30/85 6 75 200 2
22 11/26/85 6 76 200 0
23 1/12/86 6 58 200 1
Date Count Temperature Pressure Malfunction
1 4/12/81 6 66 50 0
2 11/12/81 6 70 50 1
3 3/22/82 6 69 50 0
4 11/11/82 6 68 50 0
5 4/04/83 6 67 50 0
6 6/18/82 6 72 50 0
7 8/30/83 6 73 100 0
8 11/28/83 6 70 100 0
9 2/03/84 6 57 200 1
10 4/06/84 6 63 200 1
11 8/30/84 6 70 200 1
12 10/05/84 6 78 200 0
13 11/08/84 6 67 200 0
14 1/24/85 6 53 200 2
15 4/12/85 6 67 200 0
16 4/29/85 6 75 200 0
17 6/17/85 6 70 200 0
18 7/29/85 6 81 200 0
19 8/27/85 6 76 200 0
20 10/03/85 6 79 200 0
21 10/30/85 6 75 200 2
22 11/26/85 6 76 200 0
23 1/12/86 6 58 200 1
#+end_example
Le jeu de données nous indique la date de l'essai, le nombre de joints
......@@ -85,7 +85,7 @@ data
#+end_src
#+RESULTS:
: Date Count Temperature Pressure Malfunction
: Date Count Temperature Pressure Malfunction
: 2 11/12/81 6 70 50 1
: 9 2/03/84 6 57 200 1
: 10 4/06/84 6 63 200 1
......@@ -99,7 +99,7 @@ la pression est quasiment toujours égale à 200, ce qui devrait
simplifier l'analyse.
Comment la fréquence d'échecs varie-t-elle avec la température ?
#+begin_src R :results output graphics :file "freq_temp.png" :exports both :width 600 :height 400 :session *R*
#+begin_src R :results file graphics :file "freq_temp.png" :exports both :width 600 :height 400 :session *R*
plot(data=data, Malfunction/Count ~ Temperature, ylim=c(0,1))
#+end_src
......@@ -133,10 +133,6 @@ Call:
glm(formula = Malfunction/Count ~ Temperature, family = binomial(link = "logit"),
data = data, weights = Count)
Deviance Residuals:
2 9 10 11 14 21 23
-0.3015 -0.2836 -0.2919 -0.3015 0.6891 0.6560 -0.2850
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.389528 3.195752 -0.435 0.664
......@@ -161,7 +157,7 @@ La température prévue le jour du décollage est de 31°F. Essayons
d'estimer la probabilité de dysfonctionnement des joints toriques à
cette température à partir du modèle que nous venons de construire:
#+begin_src R :results output graphics :file "proba_estimate.png" :exports both :width 600 :height 400 :session *R*
#+begin_src R :results file graphics :file "proba_estimate.png" :exports both :width 600 :height 400 :session *R*
# shuttle=shuttle[shuttle$r!=0,]
tempv = seq(from=30, to=90, by = .5)
rmv <- predict(logistic_reg,list(Temperature=tempv),type="response")
......@@ -180,7 +176,7 @@ précédents où nous il y a eu défaillance d'au moins un joint. Revenons
défaillance d'un joint:
#+begin_src R :results output :session *R* :exports both
data_full = read.csv("shuttle.csv",header=T)
data_full = read.csv("C:/Users/Jade/mooc-rr/module2/exo5/shuttle.csv",header=T)
sum(data_full$Malfunction)/sum(data_full$Count)
#+end_src
......
......@@ -70,10 +70,10 @@ Après avoir téléchargé les données, nous commençons par l'extraction des d
from urllib.request import urlopen, urlretrieve
try:
data = open('./data.csv', 'r').read()
data = open('./data.csv', 'rb').read()
except FileNotFoundError:
dataurl = urlretrieve(data_url,'./data.csv')
data = open('./data.csv', 'r').read()
data = open('./data.csv', 'rb').read()
lines = data.decode('latin-1').strip().split('\n')
data_lines = lines[1:]
......@@ -83,6 +83,7 @@ table = [line.split(',') for line in data_lines]
Regardons ce que nous avons obtenu:
#+begin_src python :results output :session :exports both
table[:5]
data.close()
#+end_src
#+RESULTS:
......
This diff is collapsed.
This diff is collapsed.
#+TITLE: Votre titre
#+AUTHOR: Votre nom
#+DATE: La date du jour
#+LANGUAGE: fr
# #+PROPERTY: header-args :eval never-export
#+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>
* Quelques explications
Ceci est un document org-mode avec quelques exemples de code
R. Une fois ouvert dans emacs, ce document peut aisément être
exporté au format HTML, PDF, et Office. Pour plus de détails sur
org-mode vous pouvez consulter https://orgmode.org/guide/.
Lorsque vous utiliserez le raccourci =C-c C-e h o=, ce document sera
compilé en html. Tout le code contenu sera ré-exécuté, les résultats
récupérés et inclus dans un document final. Si vous ne souhaitez pas
ré-exécuter tout le code à chaque fois, il vous suffit de supprimer
le # et l'espace qui sont devant le ~#+PROPERTY:~ au début de ce
document.
Comme nous vous l'avons montré dans la vidéo, on inclut du code
R de la façon suivante (et on l'exécute en faisant ~C-c C-c~):
#+begin_src R :results output :exports both
print("Hello world!")
#+end_src
#+RESULTS:
: [1] "Hello world!"
Voici la même chose, mais avec une session R (c'est le cas le
plus courant, R étant vraiment un langage interactif), donc une
persistance d'un bloc à l'autre (et on l'exécute toujours en faisant
~C-c C-c~).
#+begin_src R :results output :session *R* :exports both
summary(cars)
#+end_src
#+RESULTS:
: speed dist
: Min. : 4.0 Min. : 2.00
: 1st Qu.:12.0 1st Qu.: 26.00
: Median :15.0 Median : 36.00
: Mean :15.4 Mean : 42.98
: 3rd Qu.:19.0 3rd Qu.: 56.00
: Max. :25.0 Max. :120.00
Et enfin, voici un exemple de sortie graphique:
#+begin_src R :results output graphics :file "./cars.png" :exports results :width 600 :height 400 :session *R*
plot(cars)
#+end_src
#+RESULTS:
[[file:./cars.png]]
Vous remarquerez le paramètre ~:exports results~ qui indique que le code
ne doit pas apparaître dans la version finale du document. Nous vous
recommandons dans le cadre de ce MOOC de ne pas changer ce paramètre
(indiquer ~both~) car l'objectif est que vos analyses de données soient
parfaitement transparentes pour être reproductibles.
Attention, la figure ainsi générée n'est pas stockée dans le document
org. C'est un fichier ordinaire, ici nommé ~cars.png~. N'oubliez pas
de le committer si vous voulez que votre analyse soit lisible et
compréhensible sur GitLab.
Enfin, pour les prochains exercices, nous ne vous fournirons pas
forcément de fichier de départ, ça sera à vous de le créer, par
exemple en repartant de ce document et de le commiter vers
gitlab. N'oubliez pas que nous vous fournissons dans les ressources de
ce MOOC une configuration avec un certain nombre de raccourcis
claviers permettant de créer rapidement les blocs de code R (en
faisant ~<r~ ou ~<R~ suivi de ~Tab~).
Maintenant, à vous de jouer! Vous pouvez effacer toutes ces
informations et les remplacer par votre document computationnel.
* Exo
#+begin_src R :results output :session *R* :exports both
library(parsedate)
options(OutDec=",")
options(width=150)
#+end_src
#+RESULTS:
: Message d'avis :
: le package 'parsedate' a été compilé avec la version R 4.4.2
#+begin_src R :results output :session *R* :exports both
data = read.csv('C:/Users/Jade/mooc-rr/module3/exo2/inc-7-PAY.csv', skip=1)
head(data)
tail(data)
#+end_src
#+RESULTS:
#+begin_example
week indicator inc inc_low inc_up inc100 inc100_low inc100_up geo_insee geo_name
1 202444 7 2354 489 4219 4 1 7 FR France
2 202443 7 2130 625 3635 3 1 5 FR France
3 202442 7 2621 1246 3996 4 2 6 FR France
4 202441 7 2035 381 3689 3 1 5 FR France
5 202440 7 2125 725 3525 3 1 5 FR France
6 202439 7 2898 1333 4463 4 2 6 FR France
week indicator inc inc_low inc_up inc100 inc100_low inc100_up geo_insee geo_name
1765 199102 7 16277 11046 21508 29 20 38 FR France
1766 199101 7 15565 10271 20859 27 18 36 FR France
1767 199052 7 19375 13295 25455 34 23 45 FR France
1768 199051 7 19080 13807 24353 34 25 43 FR France
1769 199050 7 11079 6660 15498 20 12 28 FR France
1770 199049 7 1143 0 2610 2 0 5 FR France
#+end_example
#+begin_src R :results output :session *R* :exports both
na_records = apply(data, 1, function(x) any(is.na(x)))
data[na_records,]
#+end_src
#+RESULTS:
: [1] week indicator inc inc_low inc_up inc100 inc100_low inc100_up geo_insee geo_name
: <0 lignes> (ou 'row.names' de longueur nulle)
#+begin_src R :results output :session *R* :exports both
class(data$week)
class(data$inc)
#+end_src
#+RESULTS:
: [1] "integer"
: [1] "integer"
#+begin_src R :results output :session *R* :exports both
convert_week = function(w) {
ws = paste(w)
iso = paste0(substring(ws,1,4), "-W", substring(ws,5,6))
as.Date(parse_iso_8601(iso))
}
data$date = as.Date(convert_week(data$week))
class(data$date)
data = data[order(data$date),]
#+end_src
#+RESULTS:
: [1] "Date"
#+begin_src R :results output :session *R* :exports both
all(diff(data$date) == 7)
#+end_src
#+RESULTS:
: [1] TRUE
#+begin_src R :results file graphics :file (org-babel-temp-file "figure" ".png") :exports both :width 600 :height 400 :session *R*
plot(data$date, data$inc, type="l", xlab="Date", ylab="Incidence hebdomadaire")
#+end_src
#+RESULTS:
[[file:c:/Users/Jade/AppData/Local/Temp/babel-GRCje9/figurealkSMs.png]]
#+begin_src R :results file graphics :file (org-babel-temp-file "figure" ".png") :exports both :width 600 :height 400 :session *R*
with(tail(data, 200), plot(date, inc, type="l", xlab="Date", ylab="Incidence hebdomadaire"))
#+end_src
#+RESULTS:
[[file:c:/Users/Jade/AppData/Local/Temp/babel-luSDal/figureAF8kaK.png]]
#+begin_src R :results output :session *R* :exports both
pic_annuel = function(annee) {
debut = paste0(annee-1,"-09-01")
fin = paste0(annee,"-09-01")
semaines = data$date > debut & data$date <= fin
sum(data$inc[semaines], na.rm=TRUE)
}
annees <- 1992:2024
inc_annuelle = data.frame(annee = annees,
incidence = sapply(annees, pic_annuel))
head(inc_annuelle)
#+end_src
#+RESULTS:
: annee incidence
: 1 1992 834935
: 2 1993 642921
: 3 1994 662750
: 4 1995 651333
: 5 1996 564994
: 6 1997 683577
#+begin_src R :results file graphics :file (org-babel-temp-file "figure" ".png") :exports both :width 600 :height 400 :session *R*
plot(inc_annuelle, type="p", xlab="Année", ylab="Incidence annuelle")
#+end_src
#+RESULTS:
[[file:c:/Users/Jade/AppData/Local/Temp/babel-GRCje9/figureO0jGw3.png]]
#+begin_src R :results output :session *R* :exports both
head(inc_annuelle[order(-inc_annuelle$incidence),])
#+end_src
#+RESULTS:
: annee incidence
: 18 2009 841233
: 1 1992 834935
: 19 2010 834077
: 25 2016 779816
: 13 2004 778914
: 12 2003 760765
#+begin_src R :results file graphics :file (org-babel-temp-file "figure" ".png") :exports both :width 600 :height 400 :session *R*
hist(inc_annuelle$incidence, breaks=10, xlab="Incidence annuelle", ylab="Nb d'observations", main="")
#+end_src
#+RESULTS:
[[file:c:/Users/Jade/AppData/Local/Temp/babel-luSDal/figureuIwnf5.png]]
#+begin_src R :results output :session *R* :exports both
head(inc_annuelle[order(inc_annuelle$incidence),])
#+end_src
#+RESULTS:
: annee incidence
: 29 2020 221183
: 32 2023 365607
: 30 2021 377933
: 33 2024 479917
: 11 2002 515343
: 27 2018 539765
Jade@LORDI-RKHGKLV94.8436:1731010242
\ No newline at end of file
......@@ -231,3 +231,4 @@ head(inc_annuelle[order(inc_annuelle$incidence),])
: 33 2024 479917
: 11 2002 515343
: 27 2018 539765
This diff is collapsed.
\relax
\providecommand\hyper@newdestlabel[2]{}
\providecommand\HyField@AuxAddToFields[1]{}
\providecommand\HyField@AuxAddToCoFields[2]{}
\@writefile{toc}{\contentsline {section}{\numberline {1}Importation des données}{2}{section.1}\protected@file@percent }
\newlabel{sec:org71a26c3}{{1}{2}{Importation des données}{section.1}{}}
\newlabel{org7793ccf}{{1}{2}{Importation des données}{section*.2}{}}
\@writefile{toc}{\contentsline {section}{\numberline {2}Vérification de la validité des données}{3}{section.2}\protected@file@percent }
\newlabel{sec:org0b0e579}{{2}{3}{Vérification de la validité des données}{section.2}{}}
\newlabel{org1f478a7}{{2}{3}{Vérification de la validité des données}{section*.3}{}}
\newlabel{org374d44e}{{2}{3}{Vérification de la validité des données}{section*.4}{}}
\newlabel{orga886697}{{2}{3}{Vérification de la validité des données}{section*.5}{}}
\@writefile{toc}{\contentsline {section}{\numberline {3}Analyse des taux de mortalité selon les catégories fumeuses/non-fumeuses}{4}{section.3}\protected@file@percent }
\newlabel{sec:org10de194}{{3}{4}{Analyse des taux de mortalité selon les catégories fumeuses/non-fumeuses}{section.3}{}}
\@writefile{toc}{\contentsline {subsection}{\numberline {3.1}Effectifs fumeuses/non-fumeuses}{4}{subsection.3.1}\protected@file@percent }
\newlabel{sec:org2cd414c}{{3.1}{4}{Effectifs fumeuses/non-fumeuses}{subsection.3.1}{}}
\newlabel{org85e73d2}{{3.1}{4}{Effectifs fumeuses/non-fumeuses}{section*.6}{}}
\newlabel{orgaa6d737}{{3.1}{5}{Effectifs fumeuses/non-fumeuses}{section*.6}{}}
\@writefile{toc}{\contentsline {subsection}{\numberline {3.2}Effectif des femmes vivantes/mortes selon leur catégorie}{5}{subsection.3.2}\protected@file@percent }
\newlabel{sec:org5feacac}{{3.2}{5}{Effectif des femmes vivantes/mortes selon leur catégorie}{subsection.3.2}{}}
\newlabel{org4c29129}{{3.2}{5}{Effectif des femmes vivantes/mortes selon leur catégorie}{section*.7}{}}
\newlabel{orgaaf641a}{{3.2}{6}{Effectif des femmes vivantes/mortes selon leur catégorie}{section*.8}{}}
\newlabel{org82955c8}{{3.2}{7}{Effectif des femmes vivantes/mortes selon leur catégorie}{section*.8}{}}
\newlabel{org3d55922}{{3.2}{7}{Effectif des femmes vivantes/mortes selon leur catégorie}{section*.8}{}}
\@writefile{toc}{\contentsline {subsection}{\numberline {3.3}Taux de mortalité par catégorie}{8}{subsection.3.3}\protected@file@percent }
\newlabel{sec:org915f5a8}{{3.3}{8}{Taux de mortalité par catégorie}{subsection.3.3}{}}
\newlabel{org5077c05}{{3.3}{8}{Taux de mortalité par catégorie}{section*.9}{}}
\newlabel{orgbd897bf}{{3.3}{8}{Taux de mortalité par catégorie}{section*.10}{}}
\@writefile{toc}{\contentsline {section}{\numberline {4}Analyse des taux de mortalité selon les catégories fumeuses/non-fumeuses avec la notion d'âge}{9}{section.4}\protected@file@percent }
\newlabel{sec:org5389db7}{{4}{9}{Analyse des taux de mortalité selon les catégories fumeuses/non-fumeuses avec la notion d'âge}{section.4}{}}
\@writefile{toc}{\contentsline {subsection}{\numberline {4.1}Effectifs fumeuses/non-fumeuses par classe d'âge}{9}{subsection.4.1}\protected@file@percent }
\newlabel{sec:orge920994}{{4.1}{9}{Effectifs fumeuses/non-fumeuses par classe d'âge}{subsection.4.1}{}}
\newlabel{org3ee9898}{{4.1}{9}{Effectifs fumeuses/non-fumeuses par classe d'âge}{section*.11}{}}
\newlabel{org598659a}{{4.1}{10}{Effectifs fumeuses/non-fumeuses par classe d'âge}{section*.12}{}}
\newlabel{org02d534e}{{4.1}{12}{Effectifs fumeuses/non-fumeuses par classe d'âge}{section*.12}{}}
\@writefile{toc}{\contentsline {subsection}{\numberline {4.2}Effectif des femmes vivantes/mortes par classe d'âge selon leur catégorie}{12}{subsection.4.2}\protected@file@percent }
\newlabel{sec:org32ecb99}{{4.2}{12}{Effectif des femmes vivantes/mortes par classe d'âge selon leur catégorie}{subsection.4.2}{}}
\newlabel{org82404d4}{{4.2}{12}{Effectif des femmes vivantes/mortes par classe d'âge selon leur catégorie}{section*.13}{}}
\newlabel{org03099d8}{{4.2}{13}{Effectif des femmes vivantes/mortes par classe d'âge selon leur catégorie}{section*.14}{}}
\newlabel{org5e50ca7}{{4.2}{14}{Effectif des femmes vivantes/mortes par classe d'âge selon leur catégorie}{section*.15}{}}
\newlabel{org8d66551}{{4.2}{16}{Effectif des femmes vivantes/mortes par classe d'âge selon leur catégorie}{section*.16}{}}
\newlabel{org6904c61}{{4.2}{17}{Effectif des femmes vivantes/mortes par classe d'âge selon leur catégorie}{section*.16}{}}
\newlabel{org6a933ec}{{4.2}{17}{Effectif des femmes vivantes/mortes par classe d'âge selon leur catégorie}{section*.16}{}}
\newlabel{org77f3367}{{4.2}{18}{Effectif des femmes vivantes/mortes par classe d'âge selon leur catégorie}{section*.16}{}}
\newlabel{orgdd35d52}{{4.2}{18}{Effectif des femmes vivantes/mortes par classe d'âge selon leur catégorie}{section*.16}{}}
\newlabel{org50f5353}{{4.2}{19}{Effectif des femmes vivantes/mortes par classe d'âge selon leur catégorie}{section*.16}{}}
\newlabel{orge3160e5}{{4.2}{19}{Effectif des femmes vivantes/mortes par classe d'âge selon leur catégorie}{section*.16}{}}
\newlabel{orgfe0719f}{{4.2}{20}{Effectif des femmes vivantes/mortes par classe d'âge selon leur catégorie}{section*.16}{}}
\newlabel{org7025150}{{4.2}{21}{Effectif des femmes vivantes/mortes par classe d'âge selon leur catégorie}{section*.16}{}}
\newlabel{org83f3853}{{4.2}{21}{Effectif des femmes vivantes/mortes par classe d'âge selon leur catégorie}{section*.16}{}}
\newlabel{org5d6de9f}{{4.2}{22}{Effectif des femmes vivantes/mortes par classe d'âge selon leur catégorie}{section*.16}{}}
\newlabel{orgbaa6de5}{{4.2}{22}{Effectif des femmes vivantes/mortes par classe d'âge selon leur catégorie}{section*.16}{}}
\newlabel{orga6b7bf6}{{4.2}{23}{Effectif des femmes vivantes/mortes par classe d'âge selon leur catégorie}{section*.16}{}}
\@writefile{toc}{\contentsline {subsection}{\numberline {4.3}Taux de mortalité par classe d'âge selon les catégories}{23}{subsection.4.3}\protected@file@percent }
\newlabel{sec:org18578aa}{{4.3}{23}{Taux de mortalité par classe d'âge selon les catégories}{subsection.4.3}{}}
\newlabel{org124c4e0}{{4.3}{24}{Taux de mortalité par classe d'âge selon les catégories}{section*.17}{}}
\newlabel{org428872c}{{4.3}{24}{Taux de mortalité par classe d'âge selon les catégories}{section*.18}{}}
\newlabel{org4af670e}{{4.3}{25}{Taux de mortalité par classe d'âge selon les catégories}{section*.19}{}}
\newlabel{orgb4f61bf}{{4.3}{26}{Taux de mortalité par classe d'âge selon les catégories}{section*.20}{}}
\@writefile{toc}{\contentsline {section}{\numberline {5}Conclusion sur ces deux analyses}{26}{section.5}\protected@file@percent }
\newlabel{sec:orga1077b4}{{5}{26}{Conclusion sur ces deux analyses}{section.5}{}}
\gdef \@abspage@last{27}
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
\BOOKMARK [1][-]{section.1}{\376\377\000I\000m\000p\000o\000r\000t\000a\000t\000i\000o\000n\000\040\000d\000e\000s\000\040\000d\000o\000n\000n\000\351\000e\000s}{}% 1
\BOOKMARK [1][-]{section.2}{\376\377\000V\000\351\000r\000i\000f\000i\000c\000a\000t\000i\000o\000n\000\040\000d\000e\000\040\000l\000a\000\040\000v\000a\000l\000i\000d\000i\000t\000\351\000\040\000d\000e\000s\000\040\000d\000o\000n\000n\000\351\000e\000s}{}% 2
\BOOKMARK [1][-]{section.3}{\376\377\000A\000n\000a\000l\000y\000s\000e\000\040\000d\000e\000s\000\040\000t\000a\000u\000x\000\040\000d\000e\000\040\000m\000o\000r\000t\000a\000l\000i\000t\000\351\000\040\000s\000e\000l\000o\000n\000\040\000l\000e\000s\000\040\000c\000a\000t\000\351\000g\000o\000r\000i\000e\000s\000\040\000f\000u\000m\000e\000u\000s\000e\000s\000/\000n\000o\000n\000-\000f\000u\000m\000e\000u\000s\000e\000s}{}% 3
\BOOKMARK [2][-]{subsection.3.1}{\376\377\000E\000f\000f\000e\000c\000t\000i\000f\000s\000\040\000f\000u\000m\000e\000u\000s\000e\000s\000/\000n\000o\000n\000-\000f\000u\000m\000e\000u\000s\000e\000s}{section.3}% 4
\BOOKMARK [2][-]{subsection.3.2}{\376\377\000E\000f\000f\000e\000c\000t\000i\000f\000\040\000d\000e\000s\000\040\000f\000e\000m\000m\000e\000s\000\040\000v\000i\000v\000a\000n\000t\000e\000s\000/\000m\000o\000r\000t\000e\000s\000\040\000s\000e\000l\000o\000n\000\040\000l\000e\000u\000r\000\040\000c\000a\000t\000\351\000g\000o\000r\000i\000e}{section.3}% 5
\BOOKMARK [2][-]{subsection.3.3}{\376\377\000T\000a\000u\000x\000\040\000d\000e\000\040\000m\000o\000r\000t\000a\000l\000i\000t\000\351\000\040\000p\000a\000r\000\040\000c\000a\000t\000\351\000g\000o\000r\000i\000e}{section.3}% 6
\BOOKMARK [1][-]{section.4}{\376\377\000A\000n\000a\000l\000y\000s\000e\000\040\000d\000e\000s\000\040\000t\000a\000u\000x\000\040\000d\000e\000\040\000m\000o\000r\000t\000a\000l\000i\000t\000\351\000\040\000s\000e\000l\000o\000n\000\040\000l\000e\000s\000\040\000c\000a\000t\000\351\000g\000o\000r\000i\000e\000s\000\040\000f\000u\000m\000e\000u\000s\000e\000s\000/\000n\000o\000n\000-\000f\000u\000m\000e\000u\000s\000e\000s\000\040\000a\000v\000e\000c\000\040\000l\000a\000\040\000n\000o\000t\000i\000o\000n\000\040\000d\000'\000\342\000g\000e}{}% 7
\BOOKMARK [2][-]{subsection.4.1}{\376\377\000E\000f\000f\000e\000c\000t\000i\000f\000s\000\040\000f\000u\000m\000e\000u\000s\000e\000s\000/\000n\000o\000n\000-\000f\000u\000m\000e\000u\000s\000e\000s\000\040\000p\000a\000r\000\040\000c\000l\000a\000s\000s\000e\000\040\000d\000'\000\342\000g\000e}{section.4}% 8
\BOOKMARK [2][-]{subsection.4.2}{\376\377\000E\000f\000f\000e\000c\000t\000i\000f\000\040\000d\000e\000s\000\040\000f\000e\000m\000m\000e\000s\000\040\000v\000i\000v\000a\000n\000t\000e\000s\000/\000m\000o\000r\000t\000e\000s\000\040\000p\000a\000r\000\040\000c\000l\000a\000s\000s\000e\000\040\000d\000'\000\342\000g\000e\000\040\000s\000e\000l\000o\000n\000\040\000l\000e\000u\000r\000\040\000c\000a\000t\000\351\000g\000o\000r\000i\000e}{section.4}% 9
\BOOKMARK [2][-]{subsection.4.3}{\376\377\000T\000a\000u\000x\000\040\000d\000e\000\040\000m\000o\000r\000t\000a\000l\000i\000t\000\351\000\040\000p\000a\000r\000\040\000c\000l\000a\000s\000s\000e\000\040\000d\000'\000\342\000g\000e\000\040\000s\000e\000l\000o\000n\000\040\000l\000e\000s\000\040\000c\000a\000t\000\351\000g\000o\000r\000i\000e\000s}{section.4}% 10
\BOOKMARK [1][-]{section.5}{\376\377\000C\000o\000n\000c\000l\000u\000s\000i\000o\000n\000\040\000s\000u\000r\000\040\000c\000e\000s\000\040\000d\000e\000u\000x\000\040\000a\000n\000a\000l\000y\000s\000e\000s}{}% 11
This diff is collapsed.
\contentsline {section}{\numberline {1}Importation des données}{2}{section.1}%
\contentsline {section}{\numberline {2}Vérification de la validité des données}{3}{section.2}%
\contentsline {section}{\numberline {3}Analyse des taux de mortalité selon les catégories fumeuses/non-fumeuses}{4}{section.3}%
\contentsline {subsection}{\numberline {3.1}Effectifs fumeuses/non-fumeuses}{4}{subsection.3.1}%
\contentsline {subsection}{\numberline {3.2}Effectif des femmes vivantes/mortes selon leur catégorie}{5}{subsection.3.2}%
\contentsline {subsection}{\numberline {3.3}Taux de mortalité par catégorie}{8}{subsection.3.3}%
\contentsline {section}{\numberline {4}Analyse des taux de mortalité selon les catégories fumeuses/non-fumeuses avec la notion d'âge}{9}{section.4}%
\contentsline {subsection}{\numberline {4.1}Effectifs fumeuses/non-fumeuses par classe d'âge}{9}{subsection.4.1}%
\contentsline {subsection}{\numberline {4.2}Effectif des femmes vivantes/mortes par classe d'âge selon leur catégorie}{12}{subsection.4.2}%
\contentsline {subsection}{\numberline {4.3}Taux de mortalité par classe d'âge selon les catégories}{23}{subsection.4.3}%
\contentsline {section}{\numberline {5}Conclusion sur ces deux analyses}{26}{section.5}%
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