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