Initial commit

README.md

+In this project, we gather reproduction attempts from the Challenger
+study. In particular, we try to reperform some of the analysis provided
+in *Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of
+Failure* by *Siddhartha R. Dalal, Edward B. Fowlkes, Bruce Hoadley*
+published in *Journal of the American Statistical Association*, Vol. 84,
+No. 408 (Dec., 1989), pp. 945-957 and available at
+[here](https://studies2.hec.fr/jahia/webdav/site/hec/shared/sites/czellarv/acces_anonyme/OringJASA_1989.pdf)
+(here is [the official JASA
+webpage](http://www.jstor.org/stable/2290069)).
+
+On the fourth page of this article, they indicate that the maximum
+likelihood estimates of the logistic regression using only temperature
+are: $`\hat{\alpha}=5.085`$ and $`\hat{\beta}=-0.1156`$ and their
+asymptotic standard errors are $`s_{\hat{\alpha}}=3.052`$ and
+$`s_{\hat{\beta}}=0.047`$. The Goodness of fit indicated for this model
+was $`G^{2} = 18.086`$ with 21 degrees of freedom. Our goal is to
+reproduce the computation behind these values and the Figure 4 of this
+article, possibly in a nicer looking way.
+
+
+[*Here is our successful replication of Dalal et al. results using
+R*](challenger.pdf).
+
+1.  Try to **replicate the computation** from Dalal et al. In case it
+    helps, we provide you with twoimplementations of this case study
+    but we encourage you to **reimplement them by yourself** using both
+    your favourite language and an other language you do not know yet.
+    -   A [Jupyter Python3 notebook](src/Python3/challenger.ipynb)
+    -   An [Rmarkdown document](src/R/challenger.Rmd)
+2.  Then **update the following table with your own results by
+    indicating in each column:**
+	-   Language: R, Python3, Julia, Perl, C...
+	-   Language version:
+	-   Main libraries: please indicate the versions of all the loaded
+		libraries
+    -   Tool: Jupyter, Rstudio, Emacs
+	-   Operating System: Linux, Mac OS X, Windows, Android, ... along with
+      its version
+	-   $`\hat{\alpha}`$ and $`\hat{\beta}`$: Identical, Similar, Different,
+Non functional (expected values are $`5.085`$ and $`-0.1156`$)
+	-   $`s_{\hat{\alpha}}`$ and $`s_{\hat{\beta}}`$: Identical, Similar,
+			Different, Non functional (expected values are $`3.052`$ and $`0.047`$)
+	-   $`G^2`$ and degree of freedom: Identical, Similar, Different, Non
+		functional (expected values are $`18.086`$ and $`21`$).
+	-   Figure: Similar, Different, Non functional, Did not succeed
+	-   Confidence region: Identical (to [the one obtained with R](challenger.pdf)), Similar, Quite Different, Did not succeed
+ 
+| Language | Language version | Main libraries                                                | Tool    | Operating System               | $`\hat{\alpha}= 5.085`$ $`\hat{\beta} = -0.1156`$ | $`s_{\hat{\alpha}} = 3.052`$ $`s_{\hat{\beta}} = 0.047`$ | $`G^{2} = 18.086`$ $`dof = 21`$ | Figure    | Confidence Region | Link to the document                                                                                                                                                                                                                                                                                                         | Author      |
+| -------- | ---------------- | ------------------------------------------------------------- | ------- | ----------------------------   | -------------                            | ----------                                | ----------                   | --------- | ----------   | -----------------------------------------------------------------                                                                                                                                                                                                                                                                                            | ----------- |
+| R        |            3.5.1 | ggplot2 3.0.0                                                 | RStudio | Debian GNU/Linux buster/sid    | Identical                                | Identical                                 | Identical                    | Identical | Identical    | [Rmd](src/R/challenger.Rmd), [pdf](src/R/challenger_debian_alegrand.pdf)                                                                                                                                                                                                                                                                                     | A. Legrand  |
+| Python   |            3.6.4 | statsmodels 0.9.0 numpy 1.13.3 pandas 0.22.0 matplotlib 2.2.2 | Jupyter | Linux Ubuntu 4.4.0-116-generic | Identical                                | Identical                                 | Identical                    | Identical | Similar      | [ipynb](src/Python3/challenger.ipynb), [pdf](src/Python3/challenger_ubuntuMOOC_alegrand.pdf)                                                                                                                                                                                                                                                                 | A. Legrand  |
+

data/shuttle.csv

+Date,Count,Temperature,Pressure,Malfunction
+4/12/81,6,66,50,0
+11/12/81,6,70,50,1
+3/22/82,6,69,50,0
+11/11/82,6,68,50,0
+4/04/83,6,67,50,0
+6/18/82,6,72,50,0
+8/30/83,6,73,100,0
+11/28/83,6,70,100,0
+2/03/84,6,57,200,1
+4/06/84,6,63,200,1
+8/30/84,6,70,200,1
+10/05/84,6,78,200,0
+11/08/84,6,67,200,0
+1/24/85,6,53,200,2
+4/12/85,6,67,200,0
+4/29/85,6,75,200,0
+6/17/85,6,70,200,0
+7/2903/85,6,81,200,0
+8/27/85,6,76,200,0
+10/03/85,6,79,200,0
+10/30/85,6,75,200,2
+11/26/85,6,76,200,0
+1/12/86,6,58,200,1

journal.md

+# Journal
\ No newline at end of file

premier essai markdone.md

+# le titre de la premiere section
+
+## le titre de la 2eme ss section
+_italique_ *italique* 
+__gras__ **gras**
+échasse fixeé 'chasse fixe' `chasse fixe`
+--barré-- ~~barrer~~
+
+un hyperline [Élaboration et conversion de documents avec Markdown et Pandoc](https://enacit1.epfl.ch/markdown-pandoc/)
+
+<!-- un commentaire -->
+![une photo de chat](https://upload.wikimedia.org/wikipedia/commons/thumb/6/64/Collage_of_Six_Cats-02.jpg/290px-Collage_of_Six_Cats-02.jpg)
+
+liste à puce
+- +1er élément
+- 2 eme element
+- 3 eme element
+
+numérotation
+1. premier
+2. deuxieme
+3. ...
+
+Je peux imbriquer
+1. ff
+    - a
+    - b
+2. ggg
+3. on continu
+    1. au début
+    2. ensuite
+
+
+on va s'arrêter là
+
+
+
+# Nouvelle Section de Niveau 1
+
+- liste de Voitures
+    1. Porsche
+    2. Ferrari
+- Liste d'éditeurs
+    1. Notepad++
+    2. TextEdit
+
+## Nouvelle Section de Niveau 2
+
+Bonjour, c'est ma nouvelle Section de niveau 2,
+qui permet de créer un heading de niveau 2 en HTML.
+
+
+### C'est un paragraphe formé de plusieurs lignes 
+
+<p> Ligne 1 du paragraphe <br> <hr>Ligne 2 du paragraphe <br> <hr> Ligne 3 du paragraphe
+</p>
+<br>
+
+*Texte en Italique*  <br>
+
+__texte en gras__ <br>
+
+___texte en italique et gras___
+
+<br>
+[Lien vers Google](https://www.google.fr/)
+<br>
+
+![Image ](http://animtoit.fr/exemple-chien/)
+
+
+
+
+    
+
+

src/Python3/challenger_Python_org.org

+# -*- coding: utf-8 -*-
+# -*- mode: org -*-
+
+#+TITLE: Challenger - Python - Emacs - Windows 7 64 bits
+
+* Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure
+
+In this document we reperform some of the analysis provided in 
+/Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of
+Failure/ by /Siddhartha R. Dalal, Edward B. Fowlkes, Bruce Hoadley/
+published in /Journal of the American Statistical Association/, Vol. 84,
+No. 408 (Dec., 1989), pp. 945-957 and available at
+http://www.jstor.org/stable/2290069. 
+
+On the fourth page of this article, they indicate that the maximum
+likelihood estimates of the logistic regression using only temperature
+are: *$\hat{\alpha}$ = 5.085* and *$\hat{\beta}$ = -0.1156* and their
+asymptotic standard errors are *$s_{\hat{\alpha}}$ = 3.052* and
+*$s_{\hat{\beta}}$ = 0.047*. The Goodness of fit indicated for this model was
+*$G^2$ = 18.086* with *21* degrees of freedom. Our goal is to reproduce
+the computation behind these values and the Figure 4 of this article,
+possibly in a nicer looking way.
+
+** Technical information on the computer on which the analysis is run
+
+We will be using the Python 3 language using the pandas, statsmodels,
+and numpy library.
+
+#+begin_src python :results output :session :exports both
+def print_imported_modules():
+    import sys
+    for name, val in sorted(sys.modules.items()):
+        if(hasattr(val, '__version__')): 
+            print(val.__name__, val.__version__)
+#        else:
+#            print(val.__name__, "(unknown version)")
+def print_sys_info():
+    import sys
+    import platform
+    print(sys.version)
+    print(platform.uname())
+
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+import statsmodels.api as sm
+import seaborn as sns
+
+print_sys_info()
+print_imported_modules()
+#+end_src
+
+** Loading and inspecting data
+
+Let's start by reading data.
+
+#+begin_src python :results output :session :exports both
+data = pd.read_csv("https://app-learninglab.inria.fr/gitlab/moocrr-session1/moocrr-reproducibility-study/raw/master/data/shuttle.csv")
+print(data)
+#+end_src
+
+We know from our previous experience on this data set that filtering
+data is a really bad idea. We will therefore process it as such.
+
+#+begin_src python :results file :session :var matplot_lib_filename=(org-babel-temp-file "figure" ".png") :exports both
+%matplotlib inline
+pd.set_option('mode.chained_assignment',None) # this removes a useless warning from pandas
+
+data["Frequency"]=data.Malfunction/data.Count
+data.plot(x="Temperature",y="Frequency",kind="scatter",ylim=[0,1])
+plt.grid(True)
+plt.tight_layout()
+
+plt.savefig(matplot_lib_filename)
+matplot_lib_filename
+#+end_src
+
+** Logistic regression
+
+Let's assume O-rings independently fail with the same probability
+which solely depends on temperature. A logistic regression should
+allow us to estimate the influence of temperature.
+
+#+begin_src python :results output :session :exports both
+import statsmodels.api as sm
+
+data["Success"]=data.Count-data.Malfunction
+data["Intercept"]=1
+
+logmodel=sm.GLM(data['Frequency'], data[['Intercept','Temperature']], 
+                family=sm.families.Binomial(sm.families.links.logit)).fit()
+
+print(logmodel.summary())
+#+end_src
+
+The maximum likelyhood estimator of the intercept and of Temperature
+are thus *$\hat{\alpha}$ = 5.0850* and *$\hat{\beta}$ = -0.1156*. This *corresponds*
+to the values from the article of Dalal /et al./ The standard errors are 
+/$s_{\hat{\alpha}}$ = 7.477/ and /$s_{\hat{\beta}}$ = 0.115/, which is *different* from
+the *3.052* and *0.04702* reported by Dallal /et al./ The deviance is
+/3.01444/ with *21* degrees of freedom. I cannot find any value similar
+to the Goodness of fit (*$G^2$ = 18.086*) reported by Dalal /et al./ There
+seems to be something wrong. Oh I know, I haven't indicated that my
+observations are actually the result of 6 observations for each rocket
+launch. Let's indicate these weights (since the weights are always the
+same throughout all experiments, it does not change the estimates of
+the fit but it does influence the variance estimates).
+
+#+begin_src python :results output :session :exports both
+logmodel=sm.GLM(data['Frequency'], data[['Intercept','Temperature']], 
+                family=sm.families.Binomial(sm.families.links.logit),
+                var_weights=data['Count']).fit()
+
+print(logmodel.summary())
+#+end_src
+
+Good, now I have recovered the asymptotic standard errors
+*$s_{\hat{\alpha}}$ = 3.052* and *$s_{\hat{\beta}}$ = 0.047*. The Goodness of fit
+(Deviance) indicated for this model is *$G^2$ = 18.086* with *21* degrees
+of freedom (Df Residuals).
+
+*I have therefore managed to fully replicate the results of the Dalal
+/et al./ article*.
+
+** Predicting failure probability
+
+The temperature when launching the shuttle was 31°F. Let's try to
+estimate the failure probability for such temperature using our model:
+
+#+begin_src python :results output :session :exports both
+data_pred = pd.DataFrame({'Temperature': np.linspace(start=30, stop=90, num=121), 
+                          'Intercept': 1})
+data_pred['Frequency'] = logmodel.predict(data_pred)
+print(data_pred.head())
+#+end_src
+
+#+begin_src python :results file :session :var matplot_lib_filename=(org-babel-temp-file "figure" ".png") :exports both
+%matplotlib inline
+
+
+data_pred.plot(x="Temperature",y="Frequency",kind="line",ylim=[0,1])
+plt.scatter(x=data["Temperature"],y=data["Frequency"])
+plt.grid(True)
+
+plt.savefig(matplot_lib_filename)
+matplot_lib_filename
+#+end_src
+
+La fonction =logmodel.predict(data_pred)= ne fonctionne pas avec les
+dernières versions de pandas (Frequency = 1 pour toutes les températures).
+
+On peut alors utiliser le code suivant pour calculer les prédictions
+et tracer la courbe :
+
+#+begin_src python :results output :session :exports both
+# Inspiring from http://blog.yhat.com/posts/logistic-regression-and-python.html
+def logit_inv(x):
+    return(np.exp(x)/(np.exp(x)+1))
+
+data_pred['Prob']=logit_inv(data_pred['Temperature'] * logmodel.params['Temperature'] + 
+                            logmodel.params['Intercept'])
+print(data_pred.head())
+#+end_src
+
+#+begin_src python :results file :session :var matplot_lib_filename=(org-babel-temp-file "figure" ".png") :exports both
+%matplotlib inline
+data_pred.plot(x="Temperature",y="Prob",kind="line",ylim=[0,1])
+plt.scatter(x=data["Temperature"],y=data["Frequency"])
+plt.grid(True)
+
+plt.savefig(matplot_lib_filename)
+matplot_lib_filename
+#+end_src
+
+This figure is very similar to the Figure 4 of Dalal /et al./ *I have
+managed to replicate the Figure 4 of the Dalal /et al./ article.*
+
+** Computing and plotting uncertainty
+
+Following the documentation of
+[Seaborn](https://seaborn.pydata.org/generated/seaborn.regplot.html),
+I use regplot.
+
+#+begin_src python :results file :session :var matplot_lib_filename=(org-babel-temp-file "figure" ".png") :exports both
+sns.set(color_codes=True)
+plt.xlim(30,90)
+plt.ylim(0,1)
+sns.regplot(x='Temperature', y='Frequency', data=data, logistic=True)
+plt.show()
+
+plt.savefig(matplot_lib_filename)
+matplot_lib_filename
+#+end_src
+
+**I think I have managed to correctly compute and plot the uncertainty
+  of my prediction.** Although the shaded area seems very similar to
+  [the one obtained by with
+  R](https://app-learninglab.inria.fr/gitlab/moocrr-session1/moocrr-reproducibility-study/raw/5c9dbef11b4d7638b7ddf2ea71026e7bf00fcfb0/challenger.pdf),
+  I can spot a few differences (e.g., the blue point for temperature
+  63 is outside)... Could this be a numerical error ? Or a difference
+  in the statistical method ? It is not clear which one is "right".

src/R/challenger.Rmd

+---
+title: "Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure"
+author: "Arnaud Legrand"
+date: "25 October 2018"
+output: pdf_document
+---
+
+In this document we reperform some of the analysis provided in 
+*Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure* by *Siddhartha R. Dalal, Edward B. Fowlkes, Bruce Hoadley* published in *Journal of the American Statistical Association*, Vol. 84, No. 408 (Dec., 1989), pp. 945-957 and available at http://www.jstor.org/stable/2290069. 
+
+On the fourth page of this article, they indicate that the maximum likelihood estimates of the logistic regression using only temperature are: $\hat{\alpha}=5.085$ and $\hat{\beta}=-0.1156$ and their asymptotic standard errors are $s_{\hat{\alpha}}=3.052$ and $s_{\hat{\beta}}=0.047$. The Goodness of fit indicated for this model was $G^2=18.086$ with 21 degrees of freedom. Our goal is to reproduce the computation behind these values and the Figure 4 of this article, possibly in a nicer looking way.
+
+# Technical information on the computer on which the analysis is run
+We will be using the R language using the ggplot2 library.
+```{r}
+library(ggplot2)
+sessionInfo()
+```
+
+Here are the available libraries
+```{r}
+devtools::session_info()
+```
+
+
+# Loading and inspecting data
+Let's start by reading data:
+```{r}
+data = read.csv("https://app-learninglab.inria.fr/gitlab/moocrr-session1/moocrr-reproducibility-study/raw/master/data/shuttle.csv",header=T)
+data
+```
+
+We know from our previous experience on this data set that filtering data is a really bad idea. We will therefore process it as such.
+
+Let's visually inspect how temperature affects malfunction:
+```{r}
+plot(data=data, Malfunction/Count ~ Temperature, ylim=c(0,1))
+```
+
+# Logistic regression
+
+Let's assume O-rings independently fail with the same probability which solely depends on temperature. A logistic regression should allow us to estimate the influence of temperature.
+
+```{r}
+logistic_reg = glm(data=data, Malfunction/Count ~ Temperature, weights=Count, 
+                   family=binomial(link='logit'))
+summary(logistic_reg)
+```
+
+The maximum likelyhood estimator of the intercept and of Temperature are thus $\hat{\alpha}=5.0849$ and $\hat{\beta}=-0.1156$ and their standard errors are $s_{\hat{\alpha}} = 3.052$ and $s_{\hat{\beta}} = 0.04702$. The Residual deviance corresponds to the Goodness of fit $G^2=18.086$ with 21 degrees of freedom. **I have therefore managed to replicate the results of the Dalal *et al.* article**.
+
+# Predicting failure probability
+The temperature when launching the shuttle was 31°F. Let's try to 
+estimate the failure probability for such temperature using our model.:
+```{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)
+```
+
+This figure is very similar to the Figure 4 of Dalal et al. **I have managed to replicate the Figure 4 of the Dalal *et al.* article.**
+
+# Confidence on the prediction
+Let's try to plot confidence intervals with ggplot2.
+```{r, fig.height=3.3}
+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, I have a warning from ggplot2 indicating *"non-integer #successes in a binomial glm!"*. This seems fishy. Furthermore, this confidence region seems huge... It seems strange to me that the uncertainty grows so large for higher temperatures. And compared to my previous call to glm, I haven't indicated the weight which accounts for the fact that each ratio Malfunction/Count corresponds to Count observations (if someone knows how to do this...). There must be something wrong.
+
+So let's provide the "raw" data to 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)
+str(data_flat)
+```
+
+Let's check whether I obtain the same regression or not:
+```{r}
+logistic_reg_flat = glm(data=data_flat, Malfunction ~ Temperature, family=binomial(link='logit'))
+summary(logistic_reg)
+```
+Perfect. The estimates and the standard errors are the same although the Residual deviance is difference since the distance is now measured with respect to each 0/1 measurement and not to ratios. Let's use plot the regression for *data_flat* along with the ratios (*data*).
+
+```{r, fig.height=3.3}
+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()
+```
+
+This confidence interval seems much more reasonable (in accordance with the data) than the previous one. Let's check whether it corresponds to the prediction obtained when calling directly predict. Obtaining the prediction can be done directly or through the link function.
+
+Here is the "direct" (response) version I used in my very first plot:
+```{r}
+pred = predict(logistic_reg_flat,list(Temperature=30),type="response",se.fit = T)
+pred
+```
+The estimated Failure probability for 30° is thus $0.834$. However, the $se.fit$ value seems pretty hard to use as I can obviously not simply add $\pm 2 se.fit$ to $fit$ to compute a confidence interval.
+
+Here is the "link" version:
+```{r}
+pred_link = predict(logistic_reg_flat,list(Temperature=30),type="link",se.fit = T)
+pred_link
+logistic_reg$family$linkinv(pred_link$fit)
+```
+I recover $0.834$ for the estimated Failure probability at 30°. But now, going through the *linkinv* function, we can use $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))
+```
+The 95% confidence interval for our estimation is thus [0.163,0.992]. This is what ggplot2 just plotted me. This seems coherent.
+
+**I am now rather confident that I have managed to correctly compute and plot the uncertainty of my prediction.** Let's be honnest, it took me a while. My first attempts were plainly wrong (I didn't know how to do this so I trusted ggplot2, which I was misusing) and did not use the correct statistical method. I also feel confident now because this has been somehow validated by other colleagues but it will be interesting that you collect other kind of plots values that you obtained, that differ and that you would probably have kept if you didn't have a reference to compare to. Please provide us with as many versions as you can.