From 6954428ee662eea81a7d383033833b3af44a69fb Mon Sep 17 00:00:00 2001 From: Matthieu Haas Date: Tue, 9 Feb 2021 09:36:48 +0100 Subject: [PATCH] =?UTF-8?q?R=C3=A9plication=20=C3=A9tude=20Challenger=5FHA?= =?UTF-8?q?AS?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- module4/Replication__regression_haas.Rmd | 129 +++++++++++++++++++++++ 1 file changed, 129 insertions(+) create mode 100644 module4/Replication__regression_haas.Rmd diff --git a/module4/Replication__regression_haas.Rmd b/module4/Replication__regression_haas.Rmd new file mode 100644 index 0000000..bae5754 --- /dev/null +++ b/module4/Replication__regression_haas.Rmd @@ -0,0 +1,129 @@ +--- +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 -- 2.18.1