Réplication étude Challenger_HAAS

module4/Replication__regression_haas.Rmd

+---
+title: 'Réplication *"Risk Analysis of the Space Shuttle: Pre-Challenger Prediction
+  of Failure"*'
+author: "Matthieu HAAS"
+date: "05/02/2021"
+output:
+  pdf_document:
+    latex_engine: xelatex
+  html_document: default
+---
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = TRUE)
+```
+Dans ce travail, l'objectif est de répliquer certaines analyses de l'article *Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure*, disponible [ici](http://www.jstore.org/stable/2290069).
+A la quatrième page de l'article, ils présentent les meilleures estimations de la régression logistique utilisant uniquement la température comme étant : $\hat{\alpha}=5.085$ et $\hat{\beta}=−0.1156$ et leur erreur standard asymtpotique respective est de : $s_{\hat{\alpha}}=3.052$ et $s_{\hat{\beta}}= 0.047$.La qualité d'ajustement de ce modèle est de : $G^2=18.086$ avec 21 degrés de liberté. Notre objectif est de reproduire les étapes derrière ce calcul et la figure 4 de l'article, en l'améliorant si possible.
+
+## Information technique sur l'ordinateur réalisant l'analyse
+
+Nous utiliserons ici un langage R et la librairie ggplot2. 
+
+```{r}
+library(ggplot2)
+sessionInfo()
+```
+
+Voici les librairies disponibles. 
+
+```{r}
+devtools::session_info()
+```
+
+## Chargement et vérification des données
+
+La première étape est de charger les données. 
+
+```{r}
+data = read.csv("https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/raw/master/data/shuttle.csv?inline=false",header=T)
+data
+```
+
+On peut représenter l'impact de la température sur les échecs.
+
+```{r}
+plot(data=data, Malfunction/Count ~ Temperature, ylim=c(0,1))
+```
+
+
+## Régression logistique
+
+Ainsi, une régression logistique doit 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)
+```
+
+Les meilleures estimations de la régression logistique utilisant uniquement la température sont : ${\hat{\alpha}}=5.08498$ et ${\hat{\beta}}=-0.11560$ et leur erreur standard asymtpotique respective est de : $s_{\hat{\alpha}}=3.05247$ et $s_{\hat{\beta}}= 0.04702$.La qualité d'ajustement de ce modèle est de : $G^2=18.086$ avec 21 degrés de liberté. De fait, je réplique les résultats de l'article de Dalal *et al*. 
+
+## Prédiction du risque d'échec
+
+Lors du lancement, la température lors du lancement était de 31°F. Estimons la probabilité d'un échec dans ces conditions 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)
+```
+
+Cette figure est proche de la figure 4 de l'article de Dalal *et al*. Ainsi, on considère qu'elle est bien répliquée. 
+
+## Confiance de la prédiction 
+
+A l'aide de ggplot, représentons l'intervalle de confiance. 
+
+```{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()
+```
+
+On a un message d'avertissement de ggplot indiquant "non-integer #successes in a binomial glm!". C'est étrange. De plus, cet intervalle de confiance est énorme. 
+
+```{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)
+```
+
+```{r}
+logistic_reg_flat =glm(data=data_flat, Malfunction~Temperature, family=binomial(link='logit'))
+summary(logistic_reg)
+```
+
+```{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()
+```
+
+```{r}
+pred =predict(logistic_reg_flat,list(Temperature=30),type="response",se.fit = T)
+pred
+```
+
+Pour une température de 30°F, le risque d'échec est de 0.834. 
+
+```{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)
+```
+
+```{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 [0.163,0.992]. C'est ce que ggplot a représenté. 
\ No newline at end of file