Skip to content

M
mooc-rr
Source
Select Git revision
...
Target
Select Git revision
 
Commits (24)
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 6269 additions and 95 deletions
incidence-PAY-3.csv 0 → 100644
View file @ c5868a0f
This diff is collapsed.
journal/Readme.md
View file @ c5868a0f
# Journal de bord du Mooc / Mooc's logbook
# Journal de bord du Mooc *Recherche reproductible : principes méthodologiques pour une science transparente*
FR
David Pinaud, 2023-08-11
Espace réservé au journal de bord du Mooc
## Module 1
EN
Utiliser GitLab et renseigner les ReadMe avec Markdown (balisage léger) :
- Markdown Basic Syntax [link](https://www.markdownguide.org/basic-syntax/)
- [Élaboration et conversion de documents avec Markdown et Pandoc](https://www.jdbonjour.ch/cours/markdown-pandoc/) *update 2023-08-11*
Reserved for the Mooc's logbook
\ No newline at end of file
Utiliser des étiquettes et les rechercher :
- étiquettes dans des documents `Markdown` : mot-clé (encadré par des caractères comme `;`) inséré en commentaire `<!-- ;code; -->`
- rajouter des étiquettes aux métadonnées d'une image (données EXIF). [ExifTool](https://exiftool.org/)
- Moteur de recherche de bureau [DocFetcher](http://docfetcher.sourceforge.net/fr/index.html)
## Module 2 : le document computationnel
Utilisation et exercices avec soit Jupyter, RStudio ou [OrgMode](https://orgmode.org/fr/)
Exemple d'études où les résultats n'étaient pas si accessibles que ça :
- Reinhart et Rogoff 2010 *Growth in a Time of Debt*. Tableaux Exel foireux, résultats non reproductibles mais largement entérinés par les politiques
- Bennett *et al.* 2010, l'IRM et le saumon mort (algo d'imagerie à consolider)
- ...
Difficultés de mise en place :
1. manque d'informations (données, choix faits etc)
2. l'ordinateur comme source d'erreur (*point and click* des tableurs, prb de programmation, manque de rigueur et d'organisation...)
3. dimension culturelle : l'article scientifique est une simplification de la procédure (pression pour gain de temps et de place)
\ No newline at end of file
module1/exo2/fichier-markdown.md
View file @ c5868a0f
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
View file @ c5868a0f
---
title: "Votre titre"
author: "Votre nom"
date: "La date du jour"
title: "À propos du calcul de pi"
author: "David Pinaud"
date: "2023-08-16"
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 cars}
summary(cars)
```{r pir}
pi
```
Et on peut aussi aisément inclure des figures. Par exemple:
## 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 pressure, echo=FALSE}
plot(pressure)
```{r buffon}
set.seed(42)
N = 100000
x = runif(N)
theta = pi/2*runif(N)
2/(mean(x+sin(theta)>1))
```
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.
## Avec un argument “fréquentiel” de surface
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\leq 1]=\pi/4$ (voir [méthode de Monte Carlo sur Wikipedia](https://fr.wikipedia.org/wiki/Méthode_de_Monte-Carlo#Détermination_de_la_valeur_de_π)). Le code suivant illustre ce fait:
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 frequentiel}
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()
```
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 approx}
4*mean(df$Accept)
```
Maintenant, à vous de jouer! Vous pouvez effacer toutes ces informations et les remplacer par votre document computationnel.
module2/exo1/toy_document_fr.html 0 → 100644
View file @ c5868a0f
This diff is collapsed.
module2/exo2/exercice_fr.Rmd
View file @ c5868a0f
---
title: "Votre titre"
author: "Votre nom"
date: "La date du jour"
title: "Savoir faire un calcul simple soi-même"
author: "David Pianud"
date: "2023-08-18"
output: html_document
---
......@@ -10,24 +10,15 @@ output: html_document
knitr::opts_chunk$set(echo = TRUE)
```
## Quelques explications
Calculer la moyenne et l'écart-type, le min, la médiane et le max des données suivantes :
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>.
` 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)
```
Et on peut aussi aisément inclure des figures. Par exemple:
```{r pressure, echo=FALSE}
plot(pressure)
```{r calcul}
x <- 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)
mean(x)
sd(x)
min(x)
median(x)
max(x)
```
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
View file @ c5868a0f
---
title: "Votre titre"
author: "Votre nom"
date: "La date du jour"
title: "Réaliser un affichage graphique"
author: "David Pinaud"
date: 2023-08-18"
output: html_document
---
......@@ -10,24 +10,11 @@ output: html_document
knitr::opts_chunk$set(echo = TRUE)
```
## Quelques explications
## Exercice 02 (3eme partie) : Réaliser un affichage graphique
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)
```
Et on peut aussi aisément inclure des figures. Par exemple:
```{r pressure, echo=FALSE}
plot(pressure)
```{r plots, echo=FALSE}
x <- 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)
plot(x, t="l", col="blue", xlab=NA, ylab=NA)
hist(x, col="blue", xlab=NULL, ylab=NULL, main=NULL)
```
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/exo5/exo5_fr.Rmd
View file @ c5868a0f
......@@ -122,3 +122,78 @@ 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.
# Analyse par D. Pinaud, 2023-09-07
## Chargement des données
Nous commençons donc par charger ces données:
```{r}
data <- read.csv("shuttle.csv",header=T)
str(data)
```
Le jeu de données nous indique la date de l'essai, le nombre de joints
toriques mesurés (il y en a 6 sur le lançeur principal), la
température (en Farenheit) et la pression (en psi), et enfin le
nombre de dysfonctionnements relevés.
```{r}
table(data$Count)
```
Pour chaque essai, tous les joints ont été inspectés, donc pas de déséquilibre d'échantillonnage.
## Inspection graphique des données
Dysfonctionnement en fonction de la température sur l'ensemble des données :
```{r}
plot(Malfunction ~ Temperature, data)
```
On voit que le fait d'enlever les tests où il n'y a pas eu de défaillance biaise fortement la vision que l'on peut avoir de l'effet de la température.
On regarde un éventuel effet de la pression :
```{r}
plot(Malfunction ~ Pressure, data)
logistic_reg <- glm(data=data, Malfunction/Count ~ Pressure, weights=Count,
family=binomial(link='logit'))
summary(logistic_reg)
```
Pas d'évidence pour dire qu'il existe un effet de la pression.
## Estimation de l'influence de la température
Supposons que chacun des 6 joints toriques est endommagé avec la même
probabilité et indépendamment des autres et que cette probabilité ne
dépend que de la température. Si on note $p(t)$ cette probabilité, le
nombre de joints $D$ dysfonctionnant lorsque l'on effectue le vol à
température $t$ suit une loi binomiale de paramètre $n=6$ et
$p=p(t)$. Pour relier $p(t)$ à $t$, on va donc effectuer une
régression logistique.
```{r}
logistic_reg = glm(data=data, Malfunction/Count ~ Temperature, weights=Count,
family=binomial(link='logit'))
summary(logistic_reg)
```
L’estimateur le plus probable du paramètre de température est de -0.11560 (et l’erreur standard de cet estimateur est de 0.049 (SE = 0.04702), différent significativement de 0. La probabilité de dysfonctionnement des joints toriques baisse significativement quand la température augmente... __Première grosse erreur__...
## 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. comme la température du jour du lancement n'est pas comprise dans la gamme des valeurs testées, nous faisons l'hypothèse que la relation entre probabilité de défaillance et température est linéaire dans cette gamme jusqu'à 31°F (extrapolation).
```{r}
tempv <- seq(from=25, to=90, by = .5)
rmv <- predict(logistic_reg, list(Temperature=tempv), type="response", se=TRUE)
plot(tempv, rmv$fit, type="l", ylim=c(0,1), xlab="Temperature", ylab="Malfunction probability", lwd=2)
points(data=data, Malfunction/Count ~ Temperature)
lines(tempv, rmv$fit-rmv$se.fit*qt(0.975, df.residual(logistic_reg)), lty=2) # CI95% SE*qt(0.975, df)
lines(tempv, rmv$fit+rmv$se.fit*qt(0.975, df.residual(logistic_reg)), lty=2)
abline(v=31, col="red")
```
Les lignes en pointillés donnent l'intervalle de confiance de prédiction à 95%.
La probabilité de défaillance d'un joint à une température 31°F est de p=`r round(predict(logistic_reg, list(Temperature=31), type="response"), 3)`, ce qui n'est pas du tout négligeable !!!
Sachant qu'il existe un joint primaire et un joint secondaire sur chacune des trois parties du
lançeur, la probabilité de défaillance à 31°F des deux joints de la même jonction du lançeur (ce qui est suffisant pour une défaillance majeure) est de $p^2$ soit $p=$`r round((predict(logistic_reg, list(Temperature=31), type="response")^2), 3)`...
module3/exo1/analyse-syndrome-grippal.Rmd
View file @ c5868a0f
---
title: "Analyse de l'incidence du syndrôme grippal"
author: "Konrad Hinsen"
author: "David Pinaud"
output:
pdf_document:
toc: true
......@@ -43,8 +43,15 @@ 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
```{r}
data = read.csv(data_url, skip=1)
Le téléchargement se fait s'il n'existe pas une copie locale du fichier téléchargé. En effet, pour nous protéger contre une éventuelle disparition ou modification du serveur du Réseau Sentinelles, nous faisons une copie locale de ce jeux de données que nous préservons avec notre analyse. Il est inutile et même risquée de télécharger les données à chaque exécution, car dans le cas d'une panne nous pourrions remplacer nos données par un fichier défectueux. Pour cette raison, nous téléchargeons les données seulement si la copie locale n'existe pas.
```{r}
if(file.exists("incidence-PAY-3.csv"))
{
data = read.csv("incidence-PAY-3.csv")
} else {
data = read.csv(data_url, skip=1, na.strings = "-")
write.csv(data, "incidence-PAY-3.csv", row.names=FALSE)
}
```
Regardons ce que nous avons obtenu:
......
module3/exo2/exercice_fr.Rmd
View file @ c5868a0f
---
title: "Votre titre"
author: "Votre nom"
date: "La date du jour"
title: "Analyse de l'incidence de la varicelle"
author: "David Pinaud"
date: "2023-09-08"
output: html_document
---
......@@ -10,24 +10,159 @@ output: html_document
knitr::opts_chunk$set(echo = TRUE)
```
## Quelques explications
## Préparation 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>.
Les données de l'incidence de la varicelle 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:
```{r}
data_url = "http://www.sentiweb.fr/datasets/incidence-PAY-7.csv?v=4oh49"
```
Voici l'explication des colonnes donnée sur le [sur le site d'origine](https://ns.sentiweb.fr/incidence/csv-schema-v1.json):
| Nom de colonne | Libellé de colonne |
|----------------+-----------------------------------------------------------------------------------------------------------------------------------|
| `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) |
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:
La première ligne du fichier CSV est un commentaire, que nous ignorons en précisant `skip=1`.
### Téléchargement
Le téléchargement se fait s'il n'existe pas une copie locale du fichier téléchargé. En effet, pour nous protéger contre une éventuelle disparition ou modification du serveur du Réseau Sentinelles, nous faisons une copie locale de ce jeux de données que nous préservons avec notre analyse. Il est inutile et même risquée de télécharger les données à chaque exécution, car dans le cas d'une panne nous pourrions remplacer nos données par un fichier défectueux. Pour cette raison, nous téléchargeons les données seulement si la copie locale n'existe pas.
```{r}
if(file.exists("incidence-PAY-7.csv"))
{
data = read.csv("incidence-PAY-7.csv")
} else {
data = read.csv(data_url, skip=1, na.strings = "-")
write.csv(data, "incidence-PAY-7.csv", row.names=FALSE)
}
```
```{r cars}
summary(cars)
Regardons ce que nous avons obtenu:
```{r}
head(data)
tail(data)
```
Et on peut aussi aisément inclure des figures. Par exemple:
Y a-t-il des points manquants dans nos données ?
```{r}
na_records = apply(data, 1, function (x) any(is.na(x)))
data[na_records,]
```
```{r pressure, echo=FALSE}
plot(pressure)
Les deux colonnes qui nous intéressent sont `week` et `inc`. Vérifions leurs classes:
```{r}
class(data$week)
class(data$inc)
```
Ce sont des entiers, tout va bien !
### Conversion des numéros de semaine
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.
La gestion des dates est toujours un sujet délicat. Il y a un grand nombre de conventions différentes qu'il ne faut pas confondre. Notre jeux de données utilise un format que peu de logiciels savent traiter: les semaines en format [ISO-8601](https://en.wikipedia.org/wiki/ISO_8601). En `R`, il est géré par la bibliothèque [parsedate](https://cran.r-project.org/package=parsedate):
```{r}
library(parsedate)
```
Pour faciliter le traitement suivant, nous remplaçons ces semaines par les dates qui correspondent aux lundis. Voici une petite fonction qui fait la conversion pour une seule valeur:
```{r}
convert_week = function(w) {
ws = paste(w)
iso = paste0(substring(ws, 1, 4), "-W", substring(ws, 5, 6))
as.character(parse_iso_8601(iso))
}
```
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.
Nous appliquons cette fonction à tous les points, créant une nouvelle colonne `date` dans notre jeu de données:
```{r}
data$date = as.Date(convert_week(data$week))
```
Vérifions qu'elle est de classe `Date`:
```{r}
class(data$date)
```
Les points sont dans l'ordre chronologique inverse, il est donc utile de les trier:
```{r}
data = data[order(data$date),]
```
C'est l'occasion pour faire une vérification: nos dates doivent être séparées d'exactement sept jours:
```{r}
all(diff(data$date) == 7)
```
### Inspection
Regardons enfin à quoi ressemblent nos données !
```{r}
plot(data$date, data$inc, type="l", xlab="Date", ylab="Incidence hebdomadaire")
```
Maintenant, à vous de jouer! Vous pouvez effacer toutes ces informations et les remplacer par votre document computationnel.
Un zoom sur les dernières années montre mieux la localisation des pics en hiver. Le creux des incidences se trouve en été.
```{r}
with(tail(data, 200), plot(date, inc, type="l", xlab="Date", ylab="Incidence hebdomadaire"))
```
## L'incidence annuelle
### Calcul
Étant donné que le pic de l'épidémie se situe plutôt en hiver, à 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.
L'argument `na.rm=True` dans la sommation précise qu'il faut supprimer les points manquants. Ce choix est raisonnable car il n'y a qu'un seul point manquant, dont l'impact ne peut pas être très fort.
```{r}
pic_annuel = function(annee) {
debut = paste0(annee-1,"-09-01")
fin = paste0(annee,"-09-01")
semaines = data$date > debut & data$date <= fin
sum(data$inc[semaines], na.rm=TRUE)
}
```
Nous devons aussi faire attention aux premières et dernières années de notre jeux de données. Les données commencent en janvier 1991, ce qui ne permet pas de quantifier complètement le pic attribué à 1991. Nous l'enlevons donc de notre analyse. Par contre, pour une exécution en septembre 2023, les données se terminent après le 1er septembre 2023, ce qui nous permet d'inclure cette année.
```{r}
annees = 1992:2023
```
Nous créons un nouveau jeu de données pour l'incidence annuelle, en applicant la fonction `pic_annuel` à chaque année:
```{r}
inc_annuelle = data.frame(annee = annees,
incidence = sapply(annees, pic_annuel))
head(inc_annuelle)
```
### Inspection
Voici les incidences annuelles en graphique:
```{r}
plot(inc_annuelle, type="p", xlab="Année", ylab="Incidence annuelle")
```
### Identification des épidémies les plus fortes
Une liste triée par ordre décroissant d'incidence annuelle permet de plus facilement repérer les valeurs les plus élevées:
```{r}
head(inc_annuelle[order(-inc_annuelle$incidence),])
inc_annuelle[inc_annuelle$incidence==max(inc_annuelle$incidence),]
```
L'année avec l'incidence la plus faible:
```{r}
inc_annuelle[inc_annuelle$incidence==min(inc_annuelle$incidence),]
```
Enfin, un histogramme montre bien que les épidémies fortes, qui touchent environ 10% de la population française, sont assez rares: il y en eu trois au cours des 35 dernières années.
```{r}
hist(inc_annuelle$incidence, breaks=10, xlab="Incidence annuelle", ylab="Nb d'observations", main="")
```
module3/exo3/exercice_fr.Rmd
View file @ c5868a0f
---
title: "Votre titre"
author: "Votre nom"
date: "La date du jour"
output: html_document
title: 'Traitement sujet 1 : Concentration de CO~2~ dans l''atmosphère depuis 1958'
author: "David Pinaud"
date: "2023-09-12"
output:
pdf_document:
df_print: paged
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
# Objectif
Il est de présenter un document reproductible réalisé selon les modes énoncés dans le MOOC. Le sujet choisi est le n°1 _Concentration de CO~2~ dans l'atmosphère depuis 1958_ avec des mesures de la concentration de CO~2~ dans l'atmosphère effectuées à l'observatoire de Mauna Loa, Hawaii, États-Unis.
Les différents points à aborder pour l'exercice sont :
1. Réaliser un graphique qui montre une oscillation périodique superposée à une évolution systématique plus lente.
2. Séparer ces deux phénomènes, en caractérisant d'une part l'oscillation périodique et d'autre part en proposant un modèle simple de la contribution lente pour pouvoir ensuite proposer une extrapolation de cette dernière évolution jusqu'à 2025 (dans le but de pouvoir valider le modèle par des observations futures).
# Données utilisées et importation
Les données sont disponibles sur le [site internet de l'institut Scripps]( https://scrippsco2.ucsd.edu/data/atmospheric_co2/primary_mlo_co2_record.html), on télécharge le fichier
"weekly_in_situ_co2_mlo.csv" avec les observations hebdomadaires [ici](https://scrippsco2.ucsd.edu/data/atmospheric_co2/mlo.html), jusqu'à la date du 2023-09-12 (date du téléchargement). Ce fichier est en copie locale et sert maintenant de base aux analyses.
```{r download}
data_url <-
"https://scrippsco2.ucsd.edu/assets/data/atmospheric/stations/in_situ_co2/weekly/weekly_in_situ_co2_mlo.csv"
if(file.exists("weekly_in_situ_co2_mlo.csv"))
{
data <- read.csv("weekly_in_situ_co2_mlo.csv")
} else {
data <- read.csv(data_url, skip=44, header=FALSE, # on n'importe pas les 44emes premières lignes
col.names=c("WeekDate", "CO2")) # colonnes "WeekDate" et "CO2".
write.csv(data, "weekly_in_situ_co2_mlo.csv", row.names=FALSE)
}
str(data)
range(data$WeekDate)
```
Le fichier contient deux colonnes :
* `WeekDate`: la date de la semaine, ajustée à 12:00 pour le jour du milieu de chaque semaine
* `CO2`: la concentration en CO~2~ en $\mu$mole CO~2~ par mole (ppm).
# Mise en forme des données et vérification
La colonne `WeekDate` est en format `character`, nous allons la transformer en format `Date`, puis classer le `data.frame` selon cette variable par ordre croissant.
```{r}
library(parsedate)
data$WeekDate <- as.Date(data$WeekDate)
data <- data[order(data$WeekDate),]
class(data$WeekDate)
## Quelques explications
```
La colonne `WeekDate` est bien en format `Date`.
```{r}
table(is.na(data$CO2))
```
Il n'y a pas de donnée manquante dans le tableau.
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:
# Représentation graphique
On représente sur un graphique la série en entier et un extrait sur 3 ans.
```{r cars}
summary(cars)
```{r}
library(ggplot2)
library(patchwork)
p1 <- ggplot(data) +
geom_line(aes(WeekDate, CO2), col="blue") +
labs(title=bquote("Atmospheric " ~CO[2]~" concentration"),
subtitle = "Mauna Loa Observatory, Hawaii ; data from 1958 to 2023",
x="Date (weekly measures)",
y = bquote(~CO[2]~" concentration (ppm)"))
p1 <- p1 + annotate("rect", xmin=as.Date("2015-01-01"), xmax=as.Date("2017-12-31"),
ymin=396, ymax=412, fill="red", alpha=0.2)
subd <- data[data$WeekDate >= "2015-01-01" & data$WeekDate <= "2017-12-31",] # extrait sur 3 ans
p2 <- ggplot(subd) +
geom_line(aes(WeekDate, CO2), col="blue") +
geom_point(aes(WeekDate, CO2), size=0.8, fill="blue", shape=21, col="red") +
labs(x=NULL, y = NULL) +
theme(axis.text=element_text(size=5), axis.title=element_text(size=3))
p1 + inset_element(p2, 0.01, 0.5, 0.5, 0.99)
```
# Décomposition de la série
La série montre une oscillation périodique superposée à une évolution systématique plus lente. Nous allons caractériser ces deux phénomènes.
Et on peut aussi aisément inclure des figures. Par exemple:
## Tendance à long terme
On teste d'abord un ajustement linéaire (modèle le plus simple).
```{r}
lm1 <- lm(CO2 ~ as.numeric(WeekDate), data)
plot(lm1, which=1)
```
La dispersion des résidus indique une tendance non-linéaire. On ajuste alors un modèle quadratique $CO2 = a \times WeekDate^2 + b \times WeekDate + c$ qui devrait mieux coller aux données.
```{r}
lm2 <- lm(CO2 ~I(as.numeric(WeekDate)^2) + as.numeric(WeekDate), data)
plot(lm2, which=1)
knitr::kable(AIC(lm1, lm2))
```{r pressure, echo=FALSE}
plot(pressure)
```
Le modèle avec une tendance quadratique s'ajuste bien mieux aux données.
```{r}
knitr::kable(summary(lm2)$coefficients)
```
On représente alors ce modèle avec une prédiction jusqu'en milieu d'année 2025.
```{r}
newdata <- data.frame(WeekDate=seq.Date(from = as.Date("1958-03-01"),
to = as.Date("2025-07-01"), by = "weeks"))
pred <- predict(lm2, newdata = newdata, se=TRUE)
pred <- cbind(newdata, as.data.frame(pred)[,1:2])
ggplot() +
geom_line(data=data, aes(WeekDate, CO2), col="blue") +
geom_line(data=pred, aes(WeekDate, fit), col="black", linewidth=1.3) +
labs(title=bquote("Atmospheric " ~CO[2]~" concentration"),
subtitle = "Mauna Loa Observatory, Hawaii ; data from 1958 to 2023",
x="Date (weekly measures)",
y = bquote(~CO[2]~" concentration (ppm)"))
```
La concentration de CO~2~ dans l'atmosphère prédite par ce modèle pour la semaine du `r max(pred$WeekDate)` est de `r round(pred$fit[pred$WeekDate==max(pred$WeekDate)], 2)` ppm $\pm$ `r round(pred$se.fit[pred$WeekDate==max(pred$WeekDate)]*qt(0.975, df.residual(lm2)), 2)` (CI 95%).
## Variations périodiques
Nous allons travailler sur les résidus du modèle quadratique pour analyser la périodicité des variations à plus fine échelle temporelle à l'aide d'une fonction d'autocorrelation.
```{r}
data$resid.lm2 <- residuals(lm2, arg = "pearson")
acf(data$resid.lm2, lag.max=60, main = "Auto-correlation plot for residuals")
```
On observe une autocorrelation positive maximale autour du lag 52, soit 1 année (52 semaines). La périodicité de concentrations de CO~2~ dans l'atmosphère est donc d'environs 1 an.
# Conclusion
Les données de concentrations en CO~2~ dans l'atmosphère mesurées à l'observatoire de Mauna Loa (Hawaii) de 1958 à 2023 montrent des variations importantes, avec une superposition de deux phénomènes :
1. une augmentation continue positive à long terme, avec une accélaration constante. En 2025, la concentration devrait atteindre `r round(pred$fit[pred$WeekDate==max(pred$WeekDate)], 2)` ppm ;
2. une périodicité de l'ordre d'un an, probablement expliquée par des variations saisonnières dans les émissions de CO~2~.
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.
module3/exo3/weekly_in_situ_co2_mlo.csv 0 → 100644
View file @ c5868a0f
This diff is collapsed.