---
title: "regression_logistique"
author: "Meline Saubin"
date: "2 avril 2020"
output: pdf_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Analyse de risque de lancement d'une fusée

On souhaite retrouver les résultats de l'analyse de risque effectuée par [Siddhartha et al. 1989](https://www.jstor.org/stable/2290069?seq=1).

## Informations techniques

Nous utiliserons le language R :

```{r}
library(ggplot2)
sessionInfo()
```
 
 Avec les librairies disponibles :
 
```{r}
devtools::session_info()
```
 
## Jeu de données

On importe le jeu de données à l'adresse suivante : [https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/raw/master/data/shuttle.csv?inline=false](https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/raw/master/data/shuttle.csv?inline=false)

```{r}
data = read.csv("https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/raw/master/data/shuttle.csv?inline=false")
data
```

## Premiers graphiques et question scientifique

La question que l'on se pose est : y-a-t-il un lien entre la proportion de dysfonctionnements et la température, si oui lequel?
On peut donc commencer par tracer la proportion de dysfonctionnements en fonction de la température :

```{r}
# Ajout de la colonne Prop_malfunction pour la proportion de dysfonctionnements :
data$Prop_malfunction = data$Malfunction/data$Count

plot(data$Temperature, data$Prop_malfunction, xlab = "Température", 
     ylab = "Proportion de dysfonctionement")
```
On a ici l'impression que pour de faibles températures, la proportion de dysfonctionnements est plus élevée. Cette observation n'a aucune valeur statistique, il va donc falloir construire un modèle de régression logistique pour vérifier cette hypothèse.

## Régression logistique

```{r}
reg_log = glm(Prop_malfunction ~ Temperature, weights = Count, family = binomial(link="logit"), 
              data = data)

summary(reg_log)
```

Le nombre d'itérations de fisher est suffisamment faible, on peut donc observer que la température a une influence significative sur la proportion de dysfonctionnement, en choissiant un seuil de risque de 5\%.
On retrouve les informations suivantes :
$\hat{\alpha}=5.08498$ et $\hat{\beta}=-0.11560$ avec leurs erreurs standard $s_{\hat{\alpha}} = 3.05247$ et $s_{\hat{\beta}} = 0.04702$. 
La déviance résiduelle est $G^2=18.086$ avec 21 degrés de liberté.

## Courbes de prédiction

On sait que pour le lancement en question, la température était de 31°F, or ici on n'a un jeu de données que sur un intervalle de températures de 53 à 81°F, on doit donc élargir l'intervalle pour la prédiction :

```{r}
intervalle_temp = seq(from=30, to=90, by = .1)
pred = predict(reg_log,list(Temperature=intervalle_temp),type="response")
plot(intervalle_temp,pred,type="l",ylim=c(0,1), xlab = "Température", 
     ylab = "Prédiction de la probabilité de dysfonctionnement")
points(data$Prop_malfunction ~ data$Temperature)

```

On aurait donc, d'après cette prédiction, une probabilité de dysfonctionnement d'environ 0.8 à la température de 31°F. Cependant on n'a pas pris en compte l'intervalle de confiance de cette prédiction. Le code suivant s'inspire des résultats de recherches sur [un forum de discussion](https://stackoverflow.com/questions/14423325/confidence-intervals-for-predictions-from-logistic-regression) pour construire un intervalle de confiance sur une prédiction.

```{r}
intervalle_temp = seq(from=30, to=90, by = .1)

# On ajoute l'argument se.fit = TRUE pour avoir les erreurs standard, avec cette fois le type "link"
pred = predict(reg_log,list(Temperature=intervalle_temp),type="link", se.fit = TRUE)

# Pour un intervalle de confiance à 95%, 
# On suppose une distribution de loi normale des erreurs autour de la valeur prédite 
critval <- 1.96 

upr <- pred$fit + (critval * pred$se.fit)
lwr <- pred$fit - (critval * pred$se.fit)
fit <- pred$fit

# On applique ensuite la fonction de lien inversée
pred$fit <- reg_log$family$linkinv(fit)
pred$upr <- reg_log$family$linkinv(upr)
pred$lwr <- reg_log$family$linkinv(lwr)


plot(intervalle_temp,pred$fit,type="l",ylim=c(0,1), xlab = "Température", 
     ylab = "Prédiction de la probabilité de dysfonctionnement")
lines(intervalle_temp,pred$lwr,type="l",ylim=c(0,1), xlab = "Température", 
      ylab = "Prédiction de la probabilité de dysfonctionnement", col = 'red')
lines(intervalle_temp,pred$upr,type="l",ylim=c(0,1), xlab = "Température", 
      ylab = "Prédiction de la probabilité de dysfonctionnement", col = 'red')
points(data$Prop_malfunction ~ data$Temperature)

```

On  retrouve une courbe similaire à l'article, avec la proportion de dysfonctionnement, et non le nombre de dysfonctionnement. Cependant, elle n'est quand même pas identique à la courbe de l'article.