Réplication effectuée Analyse de Risque

.gitignore

+.Rproj.user
+.Rhistory
+.RData
+.Ruserdata

module4/exo5_fr_Replication.Rmd

+---
+title: "Analyse du risque de défaillance des joints toriques de la navette Challenger"
+author: "Paul Faye"
+date: "25 Mai 2021"
+output: html_document
+---
+
+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.
+
+```{r}
+library(ggplot2)
+sessionInfo()
+```
+
+
+```{r}
+devtools::session_info()
+```
+
+# Chargement des données
+
+```{r}
+data = read.csv("https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/raw/master/data/shuttle.csv",header=T)
+data
+```
+Nous savons, grâce à notre expérience précédente sur cet ensemble de données, que le filtrage des données est une très mauvaise idée.
+
+Nous allons donc les traiter comme telles, sans modification. 
+
+Vérifions visuellement comment la température affecte le dysfonctionnement :
+```{r}
+plot(data=data, Malfunction/Count ~ Temperature, ylim=c(0,1))
+```
+Supposons que les joints toriques se rompent indépendamment les uns des autres avec la même probabilité qui dépend uniquement de la température. 
+Une régression logistique devrait nous permettre d'estimer l'influence de la température.
+```{r}
+logistic_reg = glm(data=data, Malfunction/Count ~ Temperature, weights=Count,
+family=binomial(link='logit'))
+summary(logistic_reg)
+```
+L'estimateur du maximum de vraisemblance de l'ordonnée à l'origine et de la Température sont donc $\hat{\alpha}$= 5,0849 et $\hat{\beta}$ = -0,1156 et leurs erreurs type sont $s_{\hat{\alpha}}$ = 3,052 et$s_{\hat{\beta}}$  = 0,04702. La déviance résiduelle correspond à la Goodness of fit $G^{2}$ = 18,086 avec 21 degrés de liberté. 
+
+# Prévision de la probabilité de défaillance
+
+la température au moment du lancement de la navette était de 31°F. Essayons d'estimer la probabilité de défaillance pour une telle température en utilisant notre modèle :
+```{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)
+```
+Essayons de tracer les intervalles de confiance avec ggplot2.
+
+```{r}
+ggplot(data, aes(y=Malfunction/Count, x=Temperature)) + geom_point(alpha=.2, size = 2, color="blue") +
+geom_smooth(method = "glm", method.args = list(family = "binomial"), fullrange=T) +
+xlim(30,90) + ylim(0,1) + theme_bw()
+
+```
+
+Mmmh, nous avons un avertissement de ggplot2 indiquant "non-integer #successes in a binomial glm !". Cela semble louche. De plus, cette région de confiance semble énorme. . . Il semble étrange que l'incertitude augmente autant pour des températures plus élevées. Et par rapport au précédent appel à glm, nous n'avons pas indiqué le poids qui tient compte du fait que chaque rapport Malfunction/Count correspond à des observations Count (si quelqu'un sait comment faire cela. . . ). Il doit y avoir un problème.
+Fournissons donc les données "brutes" à ggplot2.
+
+```{r}
+data_flat=data.frame()
+for(i in 1:nrow(data)) {
+temperature = data[i,"Temperature"];
+malfunction = data[i,"Malfunction"];
+d = data.frame(Temperature=temperature,Malfunction=rep(0,times = data[i,"Count"]))
+if(malfunction>0) {
+d[1:malfunction, "Malfunction"]=1;
+}
+data_flat=rbind(data_flat,d)
+}
+dim(data_flat)
+```
+
+
+```{r}
+str(data_flat)
+```
+
+Vérifions si j'obtiens la même régression ou non :
+```{r}
+logistic_reg_flat = glm(data=data_flat, Malfunction ~ Temperature, family=binomial(link='logit'))
+summary(logistic_reg)
+```
+Parfait. Les estimations et les erreurs type sont les mêmes bien que la déviance résiduelle soit différente puisque la distance est maintenant mesurée par rapport à chaque mesure 0/1 et non par rapport à des ratios. Utilisons le tracé la régression pour data_flat ainsi que les ratios (données).
+
+```{r}
+ggplot(data=data_flat, aes(y=Malfunction, x=Temperature)) +
+geom_smooth(method = "glm", method.args = list(family = "binomial"), fullrange=T) +
+geom_point(data=data, aes(y=Malfunction/Count, x=Temperature),alpha=.2, size = 2, color="blue") +
+geom_point(alpha=.5, size = .5) +
+xlim(30,90) + ylim(0,1) + theme_bw()
+```
+Cet intervalle de confiance semble beaucoup plus raisonnable (en accord avec les données) que le précédent.
+Vérifions s'il correspond à la prédiction obtenue en appelant directement predict. L'obtention de la prédiction peut se faire directement ou par l'intermédiaire de la fonction link.
+Voici la version "directe" (réponse) que j'ai utilisée dans mon tout premier graphique :
+```{r}
+pred = predict(logistic_reg_flat,list(Temperature=30),type="response",se.fit = T)
+pred
+```
+La probabilité d'échec estimée pour 30° est donc de 0,834. Cependant, la valeur se.f it semble assez difficile à utiliser ainsi.
+Je ne peux évidemment pas simplement ajouter ±2 se.fit à fit pour calculer un intervalle de confiance.
+
+Voici la version "link" (lien)
+
+```{r}
+pred_link = predict(logistic_reg_flat,list(Temperature=30),type="link",se.fit = T)
+pred_link
+```
+
+```{r}
+logistic_reg$family$linkinv(pred_link$fit)
+```
+
+Je récupère 0.834 pour la probabilité d'échec estimée à 30°. Mais maintenant, en passant par la fonction linkinv, nous pouvons utiliser se.fit :
+
+```{r}
+critval = 1.96
+logistic_reg$family$linkinv(c(pred_link$fit-critval*pred_link$se.fit,
+pred_link$fit+critval*pred_link$se.fit))
+```
+
+
+L'intervalle de confiance à 95% pour notre estimation est donc de [0.163,0.992]. C'est ce que ggplot2 vient de me tracer. Cela semble cohérent.
+Je suis maintenant assez confiant dans le fait que j'ai réussi à calculer et à représenter correctement l'incertitude de ma prédiction. Soyons honnêtes, cela m'a pris du temps. Mes premières tentatives étaient manifestement erronées (je ne savais pas comment faire, alors j'ai fait confiance à ggplot2, que j'ai mal utilisé) et je n'ai pas utilisé la bonne méthode statistique. Je me sens également confiant maintenant parce que cela a été validé d'une manière ou d'une autre par d'autres collègues mais il sera intéressant que vous rassembliez d'autres types de tracés des valeurs que vous avez obtenues, qui diffèrent et que vous auriez probablement gardées si vous n'aviez pas de référence à laquelle comparer. 
+Veuillez nous fournir autant de versions que vous le pouvez.
+