journal/Journal.md

+_Jour 1_
+Ensoleillé
+
+_Jour 2_
+Nuageux
\ No newline at end of file

module1/exo2/fichier-markdown.md

+# Partie 1
+
+## Sous-partie 1 : texte
+
+Une phrase sans rien
+
+_Une phrase en italique_
+
+__Une phrase en gras__
+
+Un lien vers [fun-mooc.fr](https://www.fun-mooc.fr)
+
+Une ligne de `code`
+
+## Sous-partie 2 : listes
+
+__Liste à puces__
+
+* item
+    * sous-item
+    * sous-item
+* item
+* item
+
+__Liste numérotée__
+
+1. item
+2. item
+3. item
+
+## Sous-partie 3 : code
+
+```
+# Extrait de code
+```
\ No newline at end of file

module2/exo1/toy_document_fr.Rmd

 ---
-title: "Votre titre"
-author: "Votre nom"
-date: "La date du jour"
+title: "À propos du calcul de pi"
+author: "*Lisa Ftn*"
+date: "*07/01/2022*"
 output: html_document
 ---
 
@@ -10,24 +10,42 @@ output: html_document
 knitr::opts_chunk$set(echo = TRUE)
 ```
 
-## Quelques explications
+## En demandant à la lib maths
 
-Ceci est un document R markdown que vous pouvez aisément exporter au format HTML, PDF, et MS Word. Pour plus de détails sur R Markdown consultez <http://rmarkdown.rstudio.com>.
+Mon ordinateur m'indique que $\pi$ vaut *approximativement*
 
-Lorsque vous cliquerez sur le bouton **Knit** ce document sera compilé afin de ré-exécuter le code R et d'inclure les résultats dans un document final. Comme nous vous l'avons montré dans la vidéo, on inclue du code R de la façon suivante:
+```{r}
+pi
+```
+
+## En utilisant la méthode des aiguilles de Buffon
+
+Mais calculé avec la **méthode** des [aiguilles de Buffon](https://fr.wikipedia.org/wiki/Aiguille_de_Buffon), on obtiendrait comme **approximation** :
 
-```{r cars}
-summary(cars)
+```{r}
+set.seed(42)
+N = 100000
+x = runif(N)
+theta = pi/2*runif(N)
+2/(mean(x+sin(theta)>1))
 ```
 
-Et on peut aussi aisément inclure des figures. Par exemple:
+## Avec un argument "fréquentiel" de surface
 
-```{r pressure, echo=FALSE}
-plot(pressure)
+Sinon, une méthode plus simple à comprendre et ne faisant pas intervenir d'appel à la fonction sinus se base sur le fait que si $X \sim U(0,1)$ et $Y \sim U(0,1)$ alors $P[X^2 + Y^2 \le 1] = \pi/4$ (voir [Méthode de Monte Carlo sur Wikipedia](https://fr.wikipedia.org/wiki/M%C3%A9thode_de_Monte-Carlo#D%C3%A9termination_de_la_valeur_de_%CF%80)). Le code suivant illustre ce fait :
+
+```{r}
+set.seed(42)
+N = 1000
+df = data.frame(X = runif(N), Y = runif(N))
+df$Accept = (df$X**2 + df$Y**2 <=1)
+library(ggplot2)
+ggplot(df, aes(x=X,y=Y,color=Accept)) + geom_point(alpha=.2) + coord_fixed() + theme_bw()
 ```
 
-Vous remarquerez le paramètre `echo = FALSE` qui indique que le code ne doit pas apparaître dans la version finale du document. Nous vous recommandons dans le cadre de ce MOOC de ne pas utiliser ce paramètre car l'objectif est que vos analyses de données soient parfaitement transparentes pour être reproductibles. 
+Il est alors aisé d'obtenir une approximation (pas terrible) de $\pi$ en comptant combien de fois, en moyenne, $X^2 + Y^2$ est inférieur à 1 :
 
-Comme les résultats ne sont pas stockés dans les fichiers Rmd, pour faciliter la relecture de vos analyses par d'autres personnes, vous aurez donc intérêt à générer un HTML ou un PDF et à le commiter.
+```{r}
+4*mean(df$Accept)
+```
 
-Maintenant, à vous de jouer! Vous pouvez effacer toutes ces informations et les remplacer par votre document computationnel.

module2/exo2/exercice_fr.Rmd

 ---
-title: "Votre titre"
-author: "Votre nom"
-date: "La date du jour"
+title: "Calculs simples"
+author: "Lisa Ftn"
+date: "11/01/2022"
 output: html_document
 ---
 
@@ -10,24 +10,21 @@ output: html_document
 knitr::opts_chunk$set(echo = TRUE)
 ```
 
-## Quelques explications
+## Création du jeu de données
 
-Ceci est un document R markdown que vous pouvez aisément exporter au format HTML, PDF, et MS Word. Pour plus de détails sur R Markdown consultez <http://rmarkdown.rstudio.com>.
+On crée le vecteur contenant les données à partir des données mises à disposition dans le MOOC.
 
-Lorsque vous cliquerez sur le bouton **Knit** ce document sera compilé afin de ré-exécuter le code R et d'inclure les résultats dans un document final. Comme nous vous l'avons montré dans la vidéo, on inclue du code R de la façon suivante:
-
-```{r cars}
-summary(cars)
+```{r}
+data = c(14.0, 7.6, 11.2, 12.8, 12.5, 9.9, 14.9, 9.4, 16.9, 10.2, 14.9, 18.1, 7.3, 9.8, 10.9,12.2, 9.9, 2.9, 2.8, 15.4, 15.7, 9.7, 13.1, 13.2, 12.3, 11.7, 16.0, 12.4, 17.9, 12.2, 16.2, 18.7, 8.9, 11.9, 12.1, 14.6, 12.1, 4.7, 3.9, 16.9, 16.8, 11.3, 14.4, 15.7, 14.0, 13.6, 18.0, 13.6, 19.9, 13.7, 17.0, 20.5, 9.9, 12.5, 13.2, 16.1, 13.5, 6.3, 6.4, 17.6, 19.1, 12.8, 15.5, 16.3, 15.2, 14.6, 19.1, 14.4, 21.4, 15.1, 19.6, 21.7, 11.3, 15.0, 14.3, 16.8, 14.0, 6.8, 8.2, 19.9, 20.4, 14.6, 16.4, 18.7, 16.8, 15.8, 20.4, 15.8, 22.4, 16.2, 20.3, 23.4, 12.1, 15.5, 15.4, 18.4, 15.7, 10.2, 8.9, 21.0)
 ```
 
-Et on peut aussi aisément inclure des figures. Par exemple:
+## Calcul de la moyenne, min, max, etc.
 
-```{r pressure, echo=FALSE}
-plot(pressure)
+```{r}
+mean(data)
+max(data)
+min(data)
+median(data)
+sd(data)
 ```
 
-Vous remarquerez le paramètre `echo = FALSE` qui indique que le code ne doit pas apparaître dans la version finale du document. Nous vous recommandons dans le cadre de ce MOOC de ne pas utiliser ce paramètre car l'objectif est que vos analyses de données soient parfaitement transparentes pour être reproductibles. 
-
-Comme les résultats ne sont pas stockés dans les fichiers Rmd, pour faciliter la relecture de vos analyses par d'autres personnes, vous aurez donc intérêt à générer un HTML ou un PDF et à le commiter.
-
-Maintenant, à vous de jouer! Vous pouvez effacer toutes ces informations et les remplacer par votre document computationnel.

module2/exo3/exercice_fr.Rmd

 ---
-title: "Votre titre"
-author: "Votre nom"
-date: "La date du jour"
+title: "Calculs simples"
+author: "Lisa Ftn"
+date: "11/01/2022"
 output: html_document
 ---
 
@@ -10,24 +10,23 @@ output: html_document
 knitr::opts_chunk$set(echo = TRUE)
 ```
 
-## Quelques explications
+## Création du jeu de données
 
-Ceci est un document R markdown que vous pouvez aisément exporter au format HTML, PDF, et MS Word. Pour plus de détails sur R Markdown consultez <http://rmarkdown.rstudio.com>.
+On crée le vecteur contenant les données à partir des données mises à disposition dans le MOOC.
 
-Lorsque vous cliquerez sur le bouton **Knit** ce document sera compilé afin de ré-exécuter le code R et d'inclure les résultats dans un document final. Comme nous vous l'avons montré dans la vidéo, on inclue du code R de la façon suivante:
-
-```{r cars}
-summary(cars)
+```{r}
+data = c(14.0, 7.6, 11.2, 12.8, 12.5, 9.9, 14.9, 9.4, 16.9, 10.2, 14.9, 18.1, 7.3, 9.8, 10.9,12.2, 9.9, 2.9, 2.8, 15.4, 15.7, 9.7, 13.1, 13.2, 12.3, 11.7, 16.0, 12.4, 17.9, 12.2, 16.2, 18.7, 8.9, 11.9, 12.1, 14.6, 12.1, 4.7, 3.9, 16.9, 16.8, 11.3, 14.4, 15.7, 14.0, 13.6, 18.0, 13.6, 19.9, 13.7, 17.0, 20.5, 9.9, 12.5, 13.2, 16.1, 13.5, 6.3, 6.4, 17.6, 19.1, 12.8, 15.5, 16.3, 15.2, 14.6, 19.1, 14.4, 21.4, 15.1, 19.6, 21.7, 11.3, 15.0, 14.3, 16.8, 14.0, 6.8, 8.2, 19.9, 20.4, 14.6, 16.4, 18.7, 16.8, 15.8, 20.4, 15.8, 22.4, 16.2, 20.3, 23.4, 12.1, 15.5, 15.4, 18.4, 15.7, 10.2, 8.9, 21.0)
 ```
 
-Et on peut aussi aisément inclure des figures. Par exemple:
+## Création des graphiques demandés
 
-```{r pressure, echo=FALSE}
-plot(pressure)
+```{r}
+plot(data, type = "l", col = "blue")
 ```
+Le premier graphique est OK.
 
-Vous remarquerez le paramètre `echo = FALSE` qui indique que le code ne doit pas apparaître dans la version finale du document. Nous vous recommandons dans le cadre de ce MOOC de ne pas utiliser ce paramètre car l'objectif est que vos analyses de données soient parfaitement transparentes pour être reproductibles. 
-
-Comme les résultats ne sont pas stockés dans les fichiers Rmd, pour faciliter la relecture de vos analyses par d'autres personnes, vous aurez donc intérêt à générer un HTML ou un PDF et à le commiter.
+```{r}
+hist(data, col = "blue", xlim = c(0,25))
+```
 
-Maintenant, à vous de jouer! Vous pouvez effacer toutes ces informations et les remplacer par votre document computationnel.
+Là on obtient pas exactement le même histogramme, il faudrait probablement utiliser ggplot2...
\ No newline at end of file

module2/exo5/exo5_fr.Rmd

@@ -5,6 +5,12 @@ date: "28 juin 2018"
 output: html_document
 ---
 
+```{r load-packages, include=FALSE}
+library(ggplot2)
+library(knitr)
+```
+
+
 Le 27 Janvier 1986, veille du décollage de la navette _Challenger_, eu
 lieu une télé-conférence de trois heures entre les ingénieurs de la
 Morton Thiokol (constructeur d'un des moteurs) et de la NASA. La
@@ -93,6 +99,16 @@ plot(tempv,rmv,type="l",ylim=c(0,1))
 points(data=data, Malfunction/Count ~ Temperature)
 ```
 
+Ajout d'une courbe ggplot pour visualiser l'intervalle de confiance.
+
+```{r}
+ggplot(data, aes(x = Temperature, y = Malfunction/Count)) + geom_point(alpha = .3, size = 3) +
+  theme_bw() + 
+  geom_smooth(method = "glm", method.args = list(family = "binomial")) + xlim(50,80)
+```
+
+Deux soucis : la distrib. binomiale ne convient pas et l'intervalle de confiance est immense.
+
 Comme on pouvait s'attendre au vu des données initiales, la
 température n'a pas d'impact notable sur la probabilité d'échec des
 joints toriques. Elle sera d'environ 0.2, comme dans les essais
@@ -122,3 +138,88 @@ fiasco, l'analyse précédente comporte (au moins) un petit
 problème... Saurez-vous le trouver ? Vous êtes libre de modifier cette
 analyse et de regarder ce jeu de données sous tous les angles afin
 d'expliquer ce qui ne va pas.
+
+
+
+## Modification des données et nouvelle analyse
+
+Il n'y a aucune raison que les tests où les deux joints n'ont pas été abimés 
+soient retirés de l'étude. Pour simplifier et pouvoir appliquer un modèle logistique, 
+on passe plutôt à 0 = aucun disfonctionnement et 1 = au moins 1 disfonctionnement.
+
+```{r}
+data = read.csv("shuttle.csv",header=T)
+data[data$Malfunction>0,5] = 1
+data
+```
+
+Problème : la pression varie beaucoup si on garde tous les points de données.
+Est-ce que cela semble avoir un impact ?
+
+```{r}
+plot(data=data, Malfunction/Count ~ Pressure, ylim=c(0,1))
+```
+
+Pas vraiment concluant... Est-ce que pression et température covarient ?
+
+```{r}
+plot(data=data, Pressure ~ Temperature)
+```
+
+On va tenter de ne garder que les points avec P = 200, pour avoir un bel éventail
+de température, sans facteur confondant.
+
+```{r}
+data = data[data$Pressure == 200,]
+```
+
+
+
+Comment la fréquence d'échecs varie-t-elle avec la température ?
+```{r}
+plot(data=data, Malfunction~ Temperature, ylim=c(0,1))
+```
+
+On va pouvoir appliquer un glm pour tester, mais globalement le risque pour une température
+de 31°F semble énorme...
+
+# Estimation de l'influence de la température
+
+
+```{r}
+logistic_reg = glm(data=data, Malfunction ~ Temperature, 
+                   family=binomial(link='logit'))
+summary(logistic_reg)
+```
+
+L'estimateur le plus probable du paramètre de température est -0.2476
+et l'erreur standard de cet estimateur est de 0.122.
+
+# Estimation de la probabilité de dysfonctionnant des joints toriques
+La température prévue le jour du décollage est de 31°F. Essayons
+d'estimer la probabilité de dysfonctionnement des joints toriques à
+cette température à partir du modèle que nous venons de construire:
+
+```{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 ~ Temperature)
+```
+
+Ajout d'une courbe ggplot pour visualiser l'intervalle de confiance.
+
+```{r}
+ggplot(data, aes(x = Temperature, y = Malfunction)) + geom_point(alpha = .3, size = 3) +
+  theme_bw() + 
+  geom_smooth(method = "glm", method.args = list(family = "binomial")) + xlim(50,80)
+```
+De combien est exactement cette estimation ?
+
+```{r}
+head(rmv)
+```
+
+
+

module3/exo1/analyse-syndrome-grippal.Rmd

@@ -21,9 +21,16 @@ knitr::opts_chunk$set(echo = TRUE)
 
 ## Préparation des données
 
-Les données de l'incidence du syndrome grippal sont disponibles du site Web du [Réseau Sentinelles](http://www.sentiweb.fr/). Nous les récupérons sous forme d'un fichier en format CSV dont chaque ligne correspond à une semaine de la période demandée. Nous téléchargeons toujours le jeu de données complet, qui commence en 1984 et se termine avec une semaine récente. L'URL est:
+Les données de l'incidence du syndrome grippal sont disponibles du site Web du [Réseau Sentinelles](http://www.sentiweb.fr/). Nous les récupérons sous forme d'un fichier en format CSV dont chaque ligne correspond à une semaine de la période demandée. Nous téléchargeons toujours le jeu de données complet, qui commence en 1984 et se termine avec une semaine récente.
+Les données sont téléchargées en local avant l'analyse.
+
 ```{r}
 data_url = "http://www.sentiweb.fr/datasets/incidence-PAY-3.csv"
+data_file = "syndrome-grippal.csv"
+if (!file.exists(data_file)) {
+    download.file(data_url, data_file, method="auto")
+}
+
 ```
 
 Voici l'explication des colonnes donnée sur le [sur le site d'origine](https://ns.sentiweb.fr/incidence/csv-schema-v1.json):
@@ -42,9 +49,10 @@ Voici l'explication des colonnes donnée sur le [sur le site d'origine](https://
 | `geo_name`   | Libellé de la zone géographique (ce libellé peut être modifié sans préavis)                                                       |
 
 La première ligne du fichier CSV est un commentaire, que nous ignorons en précisant `skip=1`.
-### Téléchargement
+### Lecture du dataset (dl en local)
+
 ```{r}
-data = read.csv(data_url, skip=1)
+data = read.csv(data_file, skip=1)
 ```
 
 Regardons ce que nous avons obtenu:

module4/MonAnalyse.Rmd

+---
+title: 'Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure'
+author: "Lisa Ftn, based on Arnaud Legrand's work"
+date: "12/01/2022"
+output:
+  html_document:
+    df_print: paged
+---
+
+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 [here](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/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/raw/master/data/shuttle.csv?inline=false",header=T)
+data
+```
+
+Let's visually inspect how temperature affects malfunction:
+```{r}
+plot(data=data, Malfunction/Count ~ Temperature, ylim=c(0,1))
+```
+
+# Logistic regression
+
+I'll only work on the raw data, because I find it more elegant than working on ratios.
+Let's provide the "raw" data to R.
+```{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)
+```
+
+Computing regression parameters.
+```{r}
+logistic_reg_flat = glm(data=data_flat, Malfunction ~ Temperature, family=binomial(link='logit'))
+summary(logistic_reg_flat)
+```
+Perfect. The estimates and the standard errors are the same as in the 1989 paper.
+
+```{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(alpha=.5, size = .5) +
+  xlim(30,90) + ylim(0,1) + theme_bw()
+```
+
+Let's check whether this confidence interval 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_flat$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_flat$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.**