diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000000000000000000000000000000000000..5b6a0652566d10360493952aec6d4a4febc77083 --- /dev/null +++ b/.gitignore @@ -0,0 +1,4 @@ +.Rproj.user +.Rhistory +.RData +.Ruserdata diff --git a/module2/exo1/toy_document_fr.Rmd b/module2/exo1/toy_document_fr.Rmd index d7e628eb0bd063b0fc52fa0bcca77c01a44ad2ed..49b5ae97a24c421d77dfd89091c7e257847f193e 100644 --- a/module2/exo1/toy_document_fr.Rmd +++ b/module2/exo1/toy_document_fr.Rmd @@ -49,3 +49,5 @@ Il est alors aisé d\'obtenir une approximation (pas terrible) de $\pi$ en compt ```{r} 4*mean(df$Accept) ``` + +x=[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] diff --git a/module2/exo5/exo5_fr.html b/module2/exo5/exo5_fr.html new file mode 100644 index 0000000000000000000000000000000000000000..7ecab50601256fcc656b7cb104e4b14d1874d8ca --- /dev/null +++ b/module2/exo5/exo5_fr.html @@ -0,0 +1,557 @@ + + + + + + + + + + + + + + +Analyse du risque de défaillance des joints toriques de la navette Challenger + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + + + + + +

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:

+
data = read.csv("shuttle.csv",header=T)
+data
+
##        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
+

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 Farenheit) 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 aucun 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.

+
data = data[data$Malfunction>0,]
+data
+
##        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 ?

+
plot(data=data, Malfunction/Count ~ Temperature, ylim=c(0,1))
+

+

À 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.

+
logistic_reg = glm(data=data, Malfunction/Count ~ Temperature, weights=Count, 
+                   family=binomial(link='logit'))
+summary(logistic_reg)
+
## 
+## 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
+## 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
+

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:

+
# 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)
+

+

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:

+
data_full = read.csv("shuttle.csv",header=T)
+sum(data_full$Malfunction)/sum(data_full$Count)
+
## [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.

+
+ + + + +
+ + + + + + + + + + + + + + + diff --git a/module3/exo1/analyse-syndrome-grippal.Rmd b/module3/exo1/analyse-syndrome-grippal.Rmd index 771e78faac371f23c921f7f7aecc87f2100e9059..19895df55bcf2f4b0681fc4be1d45ed6d1d2b7d7 100644 --- a/module3/exo1/analyse-syndrome-grippal.Rmd +++ b/module3/exo1/analyse-syndrome-grippal.Rmd @@ -22,48 +22,60 @@ knitr::opts_chunk$set(echo = TRUE) ## Préparation des données Les données de l'incidence du syndrome grippal sont disponibles du site Web du [Réseau Sentinelles](http://www.sentiweb.fr/). Nous les récupérons sous forme d'un fichier en format CSV dont chaque ligne correspond à une semaine de la période demandée. Nous téléchargeons toujours le jeu de données complet, qui commence en 1984 et se termine avec une semaine récente. L'URL est: + ```{r} data_url = "http://www.sentiweb.fr/datasets/incidence-PAY-3.csv" + +# Download data from url +data_downloaded <- read.csv(data_url, skip=1) + +# Save data to csv +save(data_downloaded, file = "data_downloaded.csv") +data <- get(load("data_downloaded.csv")) ``` Voici l'explication des colonnes donnée sur le [sur le site d'origine](https://ns.sentiweb.fr/incidence/csv-schema-v1.json): | Nom de colonne | Libellé de colonne | -|----------------+-----------------------------------------------------------------------------------------------------------------------------------| -| `week` | Semaine calendaire (ISO 8601) | -| `indicator` | Code de l'indicateur de surveillance | -| `inc` | Estimation de l'incidence de consultations en nombre de cas | -| `inc_low` | Estimation de la borne inférieure de l'IC95% du nombre de cas de consultation | -| `inc_up` | Estimation de la borne supérieure de l'IC95% du nombre de cas de consultation | -| `inc100` | Estimation du taux d'incidence du nombre de cas de consultation (en cas pour 100,000 habitants) | -| `inc100_low` | Estimation de la borne inférieure de l'IC95% du taux d'incidence du nombre de cas de consultation (en cas pour 100,000 habitants) | -| `inc100_up` | Estimation de la borne supérieure de l'IC95% du taux d'incidence du nombre de cas de consultation (en cas pour 100,000 habitants) | -| `geo_insee` | Code de la zone géographique concernée (Code INSEE) http://www.insee.fr/fr/methodes/nomenclatures/cog/ | -| `geo_name` | Libellé de la zone géographique (ce libellé peut être modifié sans préavis) | - -La première ligne du fichier CSV est un commentaire, que nous ignorons en précisant `skip=1`. -### Téléchargement +|--------------|----------------------------------------------------------| +| `week` | Semaine calendaire (ISO 8601) | +| `indicator` | Code de l'indicateur de surveillance | +| `inc` | Estimation de l'incidence de consultations en nombre de cas | +| `inc_low` | Estimation de la borne inférieure de l'IC95% du nombre de cas de consultation | +| `inc_up` | Estimation de la borne supérieure de l'IC95% du nombre de cas de consultation | +| `inc100` | Estimation du taux d'incidence du nombre de cas de consultation (en cas pour 100,000 habitants) | +| `inc100_low` | Estimation de la borne inférieure de l'IC95% du taux d'incidence du nombre de cas de consultation (en cas pour 100,000 habitants) | +| `inc100_up` | Estimation de la borne supérieure de l'IC95% du taux d'incidence du nombre de cas de consultation (en cas pour 100,000 habitants) | +| `geo_insee` | Code de la zone géographique concernée (Code INSEE) | +| `geo_name` | Libellé de la zone géographique (ce libellé peut être modifié sans préavis) | + +La première ligne du fichier CSV est un commentaire, que nous ignorons en précisant `skip=1`. \### Téléchargement + ```{r} data = read.csv(data_url, skip=1) ``` Regardons ce que nous avons obtenu: + ```{r} head(data) tail(data) ``` Y a-t-il des points manquants dans nos données ? + ```{r} na_records = apply(data, 1, function (x) any(is.na(x))) data[na_records,] ``` Les deux colonnes qui nous intéressent sont `week` et `inc`. Vérifions leurs classes: + ```{r} class(data$week) class(data$inc) ``` + Ce sont des entiers, tout va bien ! ### Conversion des numéros de semaine @@ -85,21 +97,25 @@ convert_week = function(w) { ``` Nous appliquons cette fonction à tous les points, créant une nouvelle colonne `date` dans notre jeu de données: + ```{r} data$date = as.Date(convert_week(data$week)) ``` Vérifions qu'elle est de classe `Date`: + ```{r} class(data$date) ``` Les points sont dans l'ordre chronologique inverse, il est donc utile de les trier: + ```{r} data = data[order(data$date),] ``` C'est l'occasion pour faire une vérification: nos dates doivent être séparées d'exactement sept jours: + ```{r} all(diff(data$date) == 7) ``` @@ -107,11 +123,13 @@ all(diff(data$date) == 7) ### Inspection Regardons enfin à quoi ressemblent nos données ! + ```{r} plot(data$date, data$inc, type="l", xlab="Date", ylab="Incidence hebdomadaire") ``` Un zoom sur les dernières années montre mieux la localisation des pics en hiver. Le creux des incidences se trouve en été. + ```{r} with(tail(data, 200), plot(date, inc, type="l", xlab="Date", ylab="Incidence hebdomadaire")) ``` @@ -120,8 +138,8 @@ with(tail(data, 200), plot(date, inc, type="l", xlab="Date", ylab="Incidence heb ### Calcul -Étant donné que le pic de l'épidémie se situe en hiver, à cheval entre deux années civiles, nous définissons la période de référence entre deux minima de l'incidence, du 1er août de l'année $N$ au 1er août de l'année $N+1$. Nous mettons l'année $N+1$ comme étiquette sur cette année décalée, car le pic de l'épidémie est toujours au début de l'année $N+1$. Comme l'incidence de syndrome grippal est très faible en été, cette modification ne risque pas de fausser nos conclusions. -L'argument `na.rm=True` dans la sommation précise qu'il faut supprimer les points manquants. Ce choix est raisonnable car il n'y a qu'un seul point manquant, dont l'impact ne peut pas être très fort. +Étant donné que le pic de l'épidémie se situe en hiver, à cheval entre deux années civiles, nous définissons la période de référence entre deux minima de l'incidence, du 1er août de l'année $N$ au 1er août de l'année $N+1$. Nous mettons l'année $N+1$ comme étiquette sur cette année décalée, car le pic de l'épidémie est toujours au début de l'année $N+1$. Comme l'incidence de syndrome grippal est très faible en été, cette modification ne risque pas de fausser nos conclusions. L'argument `na.rm=True` dans la sommation précise qu'il faut supprimer les points manquants. Ce choix est raisonnable car il n'y a qu'un seul point manquant, dont l'impact ne peut pas être très fort. + ```{r} pic_annuel = function(annee) { debut = paste0(annee-1,"-08-01") @@ -132,11 +150,13 @@ pic_annuel = function(annee) { ``` Nous devons aussi faire attention aux premières et dernières années de notre jeux de données. Les données commencent en octobre 1984, ce qui ne permet pas de quantifier complètement le pic attribué à 1985. Nous l'enlevons donc de notre analyse. Par contre, pour une exécution en octobre 2018, les données se terminent après le 1er août 2018, ce qui nous permet d'inclure cette année. + ```{r} annees = 1986:2018 ``` Nous créons un nouveau jeu de données pour l'incidence annuelle, en applicant la fonction `pic_annuel` à chaque année: + ```{r} inc_annuelle = data.frame(annee = annees, incidence = sapply(annees, pic_annuel)) @@ -146,6 +166,7 @@ head(inc_annuelle) ### Inspection Voici les incidences annuelles en graphique: + ```{r} plot(inc_annuelle, type="p", xlab="Année", ylab="Incidence annuelle") ``` @@ -153,11 +174,13 @@ plot(inc_annuelle, type="p", xlab="Année", ylab="Incidence annuelle") ### Identification des épidémies les plus fortes Une liste triée par ordre décroissant d'incidence annuelle permet de plus facilement repérer les valeurs les plus élevées: + ```{r} head(inc_annuelle[order(-inc_annuelle$incidence),]) ``` Enfin, un histogramme montre bien que les épidémies fortes, qui touchent environ 10% de la population française, sont assez rares: il y en eu trois au cours des 35 dernières années. + ```{r} hist(inc_annuelle$incidence, breaks=10, xlab="Incidence annuelle", ylab="Nb d'observations", main="") ``` diff --git a/module3/exo1/data_downloaded.csv b/module3/exo1/data_downloaded.csv new file mode 100644 index 0000000000000000000000000000000000000000..21a50b4c5ffafeb97aa612cec74ba6fec7696e24 Binary files /dev/null and b/module3/exo1/data_downloaded.csv differ diff --git a/module3/exo1/data_url.csv b/module3/exo1/data_url.csv new file mode 100644 index 0000000000000000000000000000000000000000..a5a20b9caa3e785de39ede3ff6d0aafa964ba2da Binary files /dev/null and b/module3/exo1/data_url.csv differ diff --git a/module3/exo1/ddata_downloaded.csv b/module3/exo1/ddata_downloaded.csv new file mode 100644 index 0000000000000000000000000000000000000000..21a50b4c5ffafeb97aa612cec74ba6fec7696e24 Binary files /dev/null and b/module3/exo1/ddata_downloaded.csv differ diff --git a/mooc-rr.Rproj b/mooc-rr.Rproj new file mode 100644 index 0000000000000000000000000000000000000000..8e3c2ebc99e2e337f7d69948b93529a437590b27 --- /dev/null +++ b/mooc-rr.Rproj @@ -0,0 +1,13 @@ +Version: 1.0 + +RestoreWorkspace: Default +SaveWorkspace: Default +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: Sweave +LaTeX: pdfLaTeX