module1/exo2/fichier-markdown.md

-test
\ No newline at end of file
+# 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 à puce
+
+- item
+    - sous-item
+    - sous-item
+- item
+- item
+
+Liste numérotée
+
+1. item
+2. item
+3. item
+
+## Sous-partie 3 : code
+
+    `Extrait de code`

module2/exo1/toy_document_fr.Rmd

 ---
-title: "Votre titre"
-author: "Votre nom"
-date: "La date du jour"
+title: "À propos du calcul de pi"
+author: "Arnaud Legrand"
+date: "25 juin 2018"
 output: html_document
 ---
 
-
 ```{r setup, include=FALSE}
 knitr::opts_chunk$set(echo = TRUE)
 ```
 
-## Quelques explications
+## En demandant à la lib maths
+Mon ordinateur m'indique que $\pi$ vaut *approximativement*
 
-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>.
+```{r}
+pi
+```
 
-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:
+## 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. 
+```
 
-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.
+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 :
+```{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: "Savoir faire un calcul simple soi-même"
+author: "aaaaaaaa"
+date: "2020-06-03"
 output: html_document
 ---
 
@@ -10,24 +10,21 @@ output: html_document
 knitr::opts_chunk$set(echo = TRUE)
 ```
 
-## Quelques explications
+## Création des 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>.
-
-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)
+L'objetif de cet exercice est de calculer la moyenne et l'écart-type, le min, la médiane et le max des données suivantes :
+```{r}
+dat <- 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:
-
-```{r pressure, echo=FALSE}
-plot(pressure)
+La fonction `summary` permet d'obtenir la majorité des résultats demandés :
+```{r}
+summary(dat)
 ```
 
-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.
+On peut obtenir l'écart-type par la fonction `sd` :
+```{r}
+sd(dat)
+round(sd(dat), 2)
+```
 
-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: "Réaliser un affichage graphique"
+author: "aaaaaaaa"
+date: "2020-06-04"
 output: html_document
 ---
 
@@ -10,24 +10,40 @@ output: html_document
 knitr::opts_chunk$set(echo = TRUE)
 ```
 
-## Quelques explications
+## Création des 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>.
+L'objetif de cet exercice est de réaliser un affichage graphique (séquence plot + histogramme) des données suivantes :
+```{r}
+dat <- 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)
+```
 
-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 message=FALSE}
+library(tidyverse)
 ```
 
-Et on peut aussi aisément inclure des figures. Par exemple:
-
-```{r pressure, echo=FALSE}
-plot(pressure)
+## Sequence plot
+Les graphiques sont obtenus en utilisant la librairie `ggplot2` (inclus dans la librairie `tidyverse`).  
+Les données doivent être sous forme de dataframe pour utiliser `ggplot`.  
+```{r}
+data.frame(dat) %>% 
+  ggplot(aes(x = seq_along(dat), y = dat)) +
+  geom_line(col = "blue") +
+  scale_x_continuous(limits = c(0, 100), breaks = seq(0, 100, by = 20), expand = c(0, 0)) +
+  scale_y_continuous(limits = c(0, 25), breaks = seq(0, 25, by = 5), expand = c(0, 0)) +
+  labs(x = "", y = "") +
+  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. 
+# Histogramme
+Bien qu'utilisant le même nombre de "bins" que dans l'image attendue (10 "bins"), le résultat obtenu n'est pas identique à l'image affichée dans le MOOC.  
+```{r}
+data.frame(dat) %>% 
+  ggplot(aes(x = dat)) +
+  geom_histogram(bins = 10, fill = "blue", col = "black") +
+  scale_y_continuous(limits = c(0, 27), breaks = seq(0, 30, by = 5), expand = c(0, 0)) +
+  labs(x = "", y = "") +
+  theme_bw()
+```
 
-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.

module3/exo1/analyse-syndrome-grippal.Rmd

@@ -43,8 +43,13 @@ Voici l'explication des colonnes donnée sur le [sur le site d'origine](https://
 
 La première ligne du fichier CSV est un commentaire, que nous ignorons en précisant `skip=1`.
 ### Téléchargement
+Afin de ne pas télécharger à chaque fois le fichier de données, nous pouvons enregistrer celui-ci localement (s'il n'existe pas déjà).  
 ```{r}
-data = read.csv(data_url, skip=1)
+dest_file <- "./incidence_grippale.csv"
+if(!file.exists(dest_file)){
+  download.file(url = data_url, destfile = dest_file, method = "auto")
+}
+data = read.csv(dest_file, skip=1)
 ```
 
 Regardons ce que nous avons obtenu:

module3/exo2/exercice_fr.Rmd

 ---
-title: "Votre titre"
-author: "Votre nom"
-date: "La date du jour"
+title: "Analyse de l'incidence de la varicelle"
+author: "aaaaaaaaa"
+date: "2020-06-04"
 output: html_document
 ---
 
@@ -10,24 +10,199 @@ output: html_document
 knitr::opts_chunk$set(echo = TRUE)
 ```
 
-## Quelques explications
 
-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>.
+## Introduction
 
-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:
+Les données de l'incidence de la varicelle sont disponibles du site Web du [Réseau Sentinelles](http://www.sentiweb.fr/).  
+Celles-ci sont disponibles au format CSV. Chaque ligne du fichier correspond à une semaine de la période demandée (les premières données datent de la fin de l'année 1990, les dernières sont celles d'une semaine récente).  
 
-```{r cars}
-summary(cars)
+L'URL du fichier est la suivante :
+```{r}
+data_url <- "http://www.sentiweb.fr/datasets/incidence-PAY-7.csv" 
 ```
 
-Et on peut aussi aisément inclure des figures. Par exemple:
+Chargement des librairies utilisées dans cette étude :  
+```{r message=FALSE}
+library(tidyverse)   # Manipulation de données, graphiques, ...
+library(knitr)
+library(kableExtra)
+library(parsedate)
+```
+
+
+
+## Téléchargement des données
+
+Nous téléchargeons le fichier de données en local (si celui-ci n'existe pas). Ceci afin de nous prémunir contre un éventuel problème de connexion à ce fichier.
+```{r}
+dest_file <- "./varicelle_incidence.csv"
+if(!file.exists(dest_file)) {
+  download.file(url = data_url, destfile = dest_file, method = "auto")
+}
+```
+
+Puis nous chargeons les données depuis le fichier local (a première ligne du fichier étant un commentaire).
+```{r}
+varicelle <- read_csv(dest_file, skip = 1)
+```
+
+Les données obtenues sont de la forme :
+```{r}
+head(varicelle) %>% 
+  kable(format = "html", escape = FALSE, align ="c") %>% 
+  kable_styling(bootstrap_options = "striped", full_width = F, position = "left")
+```
+
+**Description du format du fichier :**
+
+- `week` : Semaine calendaire (ISO 8601)
+- `indicator` : Code de l'indicateur de surveillance
+- `inc` : Estimation de l'incidence de consultations en nombre de cas
+- `inc_low` : Estimation de la borne inférieure de l'IC95% du nombre de cas de consultation
+- `inc_up` : Estimation de la borne supérieure de l'IC95% du nombre de cas de consultation
+- `inc100` : Estimation du taux d'incidence du nombre de cas de consultation (en cas pour 100,000 habitants)
+- `inc100_low` : Estimation de la borne inférieure de l'IC95% du taux d'incidence du nombre de cas de consultation (en cas pour 100,000 habitants)
+- `inc100_up` : Estimation de la borne supérieure de l'IC95% du taux d'incidence du nombre de cas de consultation (en cas pour 100,000 habitants)
+- `geo_insee` : Code de la zone géographique concernée (Code INSEE) http://www.insee.fr/fr/methodes/nomenclatures/cog/
+- `geo_name` : Libellé de la zone géographique (ce libellé peut être modifié sans préavis)
+
+&nbsp;
+
+Dans un premier temps, vérifions s'il y a des données manquantes. 
+```{r}
+varicelle[rowSums(is.na(varicelle)) > 0, ]
+```
+
+Il semble que le fichier soit complet.  
+
+
+## Formattage des données
+
+### Transformation des données "date"
+
+Le format des dates est particulier dans ce fichier : il est sous une forme numérique 'aaaass' , où 'a' et 's' représentent l'année et la semaine de la mesure respectivement (selon la norme ISO 8601).  
+```{r}
+head(varicelle$week)
+tail(varicelle$week)
+```
+
+La librairie `parsedate` permet de gérer ce type de format, en modifiant légérement la valeur de la date sous la forme 'aaaa-Wss'.
+```{r}
+varicelle <- varicelle %>% 
+  # ajout du caractère "W" à la valeur de date, puis extraction en format "date"
+  mutate(temp_week = str_replace(string = week, pattern = "(\\d{4})(\\d{2})", replacement = "\\1-W\\2"),
+         iso_week = parse_iso_8601(temp_week)) %>% 
+  select(-temp_week)
+
+varicelle %>% 
+  select(week, iso_week) %>% 
+  head()
+```
+
+On peut vérifier la classe de cette nouvelle donnée :
+```{r}
+class(varicelle$iso_week)
+```
+
+Enfin, on peut trier les données afin de les avoir en ordre chronologique.  
+```{r}
+varicelle <- varicelle %>% 
+  arrange(iso_week)
+```
+
+### Visualisation
 
-```{r pressure, echo=FALSE}
-plot(pressure)
+Une première analyse de l'incidence de la varicelle peut être faite visuellement.  
+```{r}
+ggplot(varicelle, aes(x = iso_week, y = inc)) +
+  geom_line(col = "blue", alpha = 0.6) +
+  labs(title = "Incidence (hebdomadaire) de la varicelle",
+       x = "Semaine", y = "Nb de cas par semaine") +
+  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 semble y avoir une certaine saisonnalité. On peut zoomer sur les dernières années pour avoir un meilleur aperçu.  
+```{r}
+varicelle %>% 
+  filter(row_number() >= length(week) - 200) %>% 
+  ggplot(aes(x = iso_week, y = inc)) +
+  geom_line(col = "blue", alpha = 0.6) +
+  labs(title = "Incidence (hebdomadaire) de la varicelle",
+       x = "Semaine", y = "Nb de cas par semaine") +
+  scale_x_datetime(date_breaks = "6 months", date_labels = "%Y-%m") +
+  theme_bw()
+```
+
+Les pics d'incidence semblent avoir lieu au mois d'avril chaque année, alors que les creux sont plutôt présents vers le mois d'août.  
+
 
-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.
+## Analyse de l'incidence
+
+### Modification de la période de référence
+
+Étant donné que le pic de l'épidémie semble se situer au mois d'avril, à cheval entre deux années civiles, nous définissons la période de référence entre deux minima de l'incidence, du 1er septembre de l'année $N$ au 1er septembre de l'année $N+1$.  
+Nous mettons l'année $N+1$ comme étiquette sur cette année décalée, car le pic de l'épidémie est toujours au début de l'année $N+1$. Comme l'incidence de la varicelle est très faible en été, cette modification ne risque pas de fausser nos conclusions. 
+```{r}
+varicelle <- varicelle %>% 
+  # extraire l'année "civile" pour chaque ligne
+  mutate(annee_civile = as.numeric(format(iso_week, format = "%Y"))) %>% 
+  group_by(annee_civile) %>% 
+  # si la date est avant le 1er septembre, année précédente, sinon année en cours
+  mutate(annee_ref = ifelse(iso_week < as.POSIXct(paste0(annee_civile, "-09-01"), format = "%Y-%m-%d"),
+                            annee_civile, annee_civile + 1)) %>% 
+  ungroup()
+```
+
+On peut regarder le résultat en extrayant quelques données autour du 1er septembre 2005.  
+```{r}
+varicelle %>% 
+  filter(between(iso_week, as.POSIXct("2005-08-10"), as.POSIXct("2005-09-20"))) %>% 
+  select(week, iso_week, annee_civile, annee_ref)
+```
+
+Enfin, l'année en cours n'étant pas complète, on peut supprimer les données liées à celle-ci.  
+```{r}
+varicelle <- varicelle %>% 
+  filter(annee_ref != max(annee_ref))
+```
+
+
+### Visualisation
+
+On peut maintenant observer l'incidence par année de référence.  
+```{r}
+var_tableau <- varicelle %>%
+  group_by(annee_ref) %>%
+  summarise(inc_totale = sum(inc)) %>%
+  ungroup() %>% 
+  mutate(inc_totale_format = cell_spec(inc_totale, 
+                                       format = "html",
+                                       bold = TRUE,
+                                       color = "white",
+                                       background = ifelse(inc_totale == min(inc_totale), "lightgreen",
+                                                           ifelse(inc_totale == max(inc_totale), "red", "white"))))
+
+var_tableau %>%
+  filter(inc_totale %in% c(min(inc_totale), max(inc_totale))) %>% 
+  select(`Année` = annee_ref, `Incidence totale` = inc_totale_format) %>% 
+  kable(format = "html", escape = FALSE, align ="c") %>%
+  kable_styling(bootstrap_options = c("striped", "bordered"), full_width = F, position = "left")
+```
+
+On retrouve les résultats du tableau dans le graphe suivant :
+```{r}
+ggplot(var_tableau, aes(x = annee_ref, y = inc_totale)) +
+  geom_line(col = "blue", alpha = 0.6) +
+  geom_label(data = var_tableau %>% filter(inc_totale == min(inc_totale)),
+             aes(x = annee_ref, y = inc_totale, label = paste("Année", annee_ref, ":", inc_totale, "cas")),
+             color = "green",
+             nudge_y = -20000) +
+  geom_label(data = var_tableau %>% filter(inc_totale == max(inc_totale)),
+             aes(x = annee_ref, y = inc_totale, label = paste("Année", annee_ref, ":", inc_totale, "cas")),
+             color = "red",
+             nudge_y = 20000) +
+  labs(title = "Incidence (annuelle) de la varicelle",
+       x = "Année de référence", y = "Nb de cas par année") +
+  theme_bw()
+```
 
-Maintenant, à vous de jouer! Vous pouvez effacer toutes ces informations et les remplacer par votre document computationnel.

module3/exo3/wheat.csv

+"","Year","Wheat","Wages"
+"1",1565,41,5
+"2",1570,45,5.05
+"3",1575,42,5.08
+"4",1580,49,5.12
+"5",1585,41.5,5.15
+"6",1590,47,5.25
+"7",1595,64,5.54
+"8",1600,27,5.61
+"9",1605,33,5.69
+"10",1610,32,5.78
+"11",1615,33,5.94
+"12",1620,35,6.01
+"13",1625,33,6.12
+"14",1630,45,6.22
+"15",1635,33,6.3
+"16",1640,39,6.37
+"17",1645,53,6.45
+"18",1650,42,6.5
+"19",1655,40.5,6.6
+"20",1660,46.5,6.75
+"21",1665,32,6.8
+"22",1670,37,6.9
+"23",1675,43,7
+"24",1680,35,7.3
+"25",1685,27,7.6
+"26",1690,40,8
+"27",1695,50,8.5
+"28",1700,30,9
+"29",1705,32,10
+"30",1710,44,11
+"31",1715,33,11.75
+"32",1720,29,12.5
+"33",1725,39,13
+"34",1730,26,13.3
+"35",1735,32,13.6
+"36",1740,27,14
+"37",1745,27.5,14.5
+"38",1750,31,15
+"39",1755,35.5,15.7
+"40",1760,31,16.5
+"41",1765,43,17.6
+"42",1770,47,18.5
+"43",1775,44,19.5
+"44",1780,46,21
+"45",1785,42,23
+"46",1790,47.5,25.5
+"47",1795,76,27.5
+"48",1800,79,28.5
+"49",1805,81,29.5
+"50",1810,99,30
+"51",1815,78,NA
+"52",1820,54,NA
+"53",1821,54,NA

module3/exo3/wheat_monarchs.csv

+"name","start","end","commonwealth"
+"Elizabeth","1565","1603",""
+"James I","1603","1625",""
+"Charles I","1625","1649",""
+"Cromwell","1649","1660","True"
+"Charles II","1660","1685",""
+"James II","1685","1689",""
+"W&M","1689","1702",""
+"Anne","1702","1714",""
+"George I","1714","1727",""
+"George II","1727","1760",""
+"George III","1760","1820",""
+"George IV","1820","1821",""

module4/challenger.Rmd

+---
+title: "Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure"
+author: "Tegegne"
+date: "June 16, 2020"
+output: html_document
+---
+
+```{r}
+library(ggplot2)
+sessionInfo()
+```
+
+Here are the available libraries
+```{r}
+devtools::session_info()
+```
+
+
+# Loading and inspecting data
+```{r}
+data <- read.csv("~/R/data.csv")
+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.