Le niveau d’etiage de la fontaine des Capucins

Catacombes de Paris

Introduction à R

Cette page a pour but de faire des petits tutos en francais pour faire des stat facilement avec R.

Les données seront cataphiles, parceque sinon ce n’est pas drôle. Tous types de données, dans un premier temps, la hauteur d’étiage de la fontaine des Capucins, puis dans un deuxième, la profondeur de puits de carrières. Si vous avez d’autres idées / données non confidentielles, je prends.

R est un logiciel absolument topissime, libre (GNU General public license, Free as in freedom), pour faire des statistiques, jouer avec des données et représenter graphiquement des données, y compris de la cartographie. Ca peut servir… Ca tourne sur les pleins de plateformes, (ex: Windows, Linux et MacOS). On peut tous faire de la petite analyse sans prétention pour s’amuser, pour un rapport de stage, ou pour des gros calculs, analyser et décrypter le génome…

Tout est paramétrable, les possibilités sont infinies.
Il existe une grosse communauté d’utilisateurs qui permet d’avoir de la documentation, des conseils et de nombreux forum actifs, et une grosse communauté de contributeurs, qui développent des packages pour des usages particuliers. Il ne faut pas hésiter à copier coller une erreur de la console dans un moteur de recherche, bien souvent d’autres ont déjà eu exactement les même bugs que vous.

Qu’est ce que la statistique ?
L’ensemble des principes et des méthodes mis au point pour recueillir, classer, synthéthiser et communiquer des données numériques en vue d’une utilisation de celles-ci
Sanders et al., 1984

A quoi servent les statistiques ?

Les stat ne disent pas la vérité.
Les stat ne compensent pas une lacune de données.
Les stat sont toujours associées à une probabilité d’erreur, et permettent de quantifier les incertitudes.

Les stat aident à obtenir des critères objectifs pour tester des hypothèses (pas de choix subjectifs par des cerveaux humains)
Les stat servent à décrire l’inconnu.

Les stat permettent d’optimiser des efforts.

On connait un échantillon, que déduire de la population dont il est extrait ?

Installation R

Pour l’installer, télécharger la dernière version ici :
https://www.r-project.org/

Puis dans un 2eme temps installer l’interface Rstudio :
https://rstudio.com/products/rstudio/

C’est pas obligatoire, mais c’est moins austère que R tout seul, il y a quelques menus contextuels, et permet le Rmarkdown (les résultats sont des jolies pages html comme celle-là) … Oui parceque R c’est un langage de programmation. Mais pas de panique. Ca veut juste dire qu’il y a un script, on écrit dessus des commandes (ou on copie colle celles que je propose là), et on les exécute (Ctrl R sur la ligne ciblée, par exemple). Si c’est un script R le résultat s’affiche dans la console, ou dans la fenêtre plot. Si c’est un script Rmardown (Rmd) on peut obtenir le résultat de tout le script dans un fichier html (cliquer sur “knit to html”), pdf, latex, pwpt, word, ou autres…

Donc à cette étape R et Rstudio sont installés.

Maintenant ouvrir Rstudio et créer un fichier. File/New file/Rmarkdown

R R R

Parametres d’affichage knit

Dans l’entete document vous avez encadré par trois tirets ceci :
---    
title: "R Notebook"    
output: html_notebook    
---    

Je vous conseille de copier coller à la place de ca, bêtement mon paramétrage, et de personnaliser à votre sauce plus tard.

Ouvrir mon script Rmd d’origine #ici

Et copié ce qui est entre les 2 x trois tirets au début (en changeant juste les title et author)

extrait :
---    
title: "Kta_Stat_R"  
author: "Kafka"  
[...]  
editor_options:   
chunk_output_type: console     
transition: linear   
---   

Sommaire cliquable et textes

Il vous suffit de mettre des 1, 2 dieses (# ou ##) pour des titres et 3 pour sous-titre (###), 4 pour les sous-sous titres (####) etc. Les titres reconnus s’affichent en bleu.

Vous pouvez écrire du texte entre chaque titre.
Pour un retour à la ligne entre deux phrases ajouter 2 espaces à la fin de la ligne.

Chunks R

Voila, maintenant écrivons une commandes R dans un chunks R.

Un chunk est un morceau de code inséré dans la page de texte, dont on va pouvoir exécuter les commandes.

Pour en ajouter un, cliquez sur le bouton vert en haut à droite insert / R chunk

chunkadd

chunkadd

Dans la page finale, on peut en afficher juste le code, uniquement le résultat (output) du code exécuté, ou les 2, ou aucun (Show nothing). Choisir avec la petite icône paramètre (la petite roue d’engrenage en haut a droite du chunk, Modify chunk option).

Pour le moment laissez par défaut. Mais cela peut etre utile de cocher “show nothing (don’t run code)” lorsque un chunk comporte une erreur et que vous voulez vérifier si le reste fonctionne.

Un autre moyen pour qu’une ligne du chunk en particulier ne soit pas exécuté on peut la commenter (en mettant un # devant). Pour commenter ou dé-commenter une ligne le raccourci clavier est Ctrl+Maj+C.

Installation packages

copiez le contenu du chunk suivant :

install.packages("ggcorrplot")
install.packages("magrittr")
install.packages("ggpubr")
install.packages("ggcorrplot")
install.packages("factoextra")
install.packages("kableExtra")

Lors de la premiere utilisation d’un package vous devez faire install.packages(nom.du.package), puis le charger avec library(nom.du.package).

A la seconde utilisation, library(nom.du.package) suffit à le charger.

Charger les packages

# charger les librairies utiles


library("knitr")

knitr::opts_chunk$set(
  fig.width = 5, fig.height = 5, 
  fig.path = 'figures/StressNet_',
  fig.align = "center", 
  size = "tiny", 
  echo = TRUE, eval = TRUE, 
  warning = FALSE, message = FALSE, 
  results = TRUE, comment = "")

library("ggplot2")
library("magrittr")
library("ggpubr")
library("ggcorrplot")
library("factoextra")
library("kableExtra")

Pour tester l’exécution du chunk cliquez sur le petit triangle vert de lecture en haut à droite de celui ci et voir le résultat dans la console.

Voila les librairies utiles pour ce que je veux faire apres sont chargées. Il existe des libraries pour toutes types d’applications. Ici c’est juste pour de la représentation graphique de donnée (ggplot2), et des représentations de correlations (ggcorrplot). magrittr et ggpubr permettent de mettre facilement les test de significativité sur graphiques type boxplot.

Maintenant il faut définir le dossier de travail, celui où vous allez mettre des fichiers de données que vous aller utiliser, donc dans le chunk ci-dessous je met mon chemin, vous devez mettre le votre.

dir.base <- "D:/Users/mferrand/Documents/R/kta"

Charger le fichier de données

Le données peuvent etre sous different formats pour etre importés, ici j’importe des fichiers excel au format csv (Comma-separated values)

Important ! les valeurs décimales doivent comporter un point “.” et pas une virgule “,”
et éviter les accents, les virgules dans les entêtes de colonnes ou de lignes, et les espaces + unités dans les cases avec des valeurs. Evitez aussi les cases vides. Pour les données manquantes écrire ‘na’. mettez bien une entete de colonne pour chaque colonne.

Ici, je charge un fichier de hauteur d’eau à l’échelle d’étiage de la fontaine des Capucins (données provenant une liaison seadacc de 2015), et des données météos (trouvées sur l’historique de méteofrance.fr).

Comme c’est sympa d’être sur R Mark down on peut même ajouter des images !

fontaine_des_capucins
Source : FlickR Hugo Clément
Ici le code pour ajouter une image dans le texte à mettre direct dans le texte (pas dans un chunk) :
![legende](http://liens.com/fichier.jpg){width=300px}

Si vous voulez récupérer mon fichier source .Csv pour tester, le voila : #ici

Je nomme ce fichier importé File1. La flèche “<-” signifit qu’on stock un objet dans une variable. ici j’ai mis header=TRUE pour signifié qu’il y a bien des noms de colonne dans la premiere ligne. Je n’ai pas de nom de ligne, mais sinon il aurait fallut ajouter row.names = TRUE. Pour les fichiers csv les séparateurs entre les cases sont des points-virgule, je les désigne par sep=“;”

pathFile1 <- file.path(dir.base, 'accu_precip_nappe.csv')

File1 <- read.table(file = pathFile1, header = TRUE, as.is = TRUE, sep = ";")

La commande head(“nom_tableau”) permet de voir un aperçu des premieres lignes du tableau. Mais pour la sortie Rmarkd c’est plus joli de faire un : kable(head(File1))

## la version basique
# kable(head(File1))

# la version parametrable
kable(head(File1)) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
date etiage_m temp_mois_min temp_mois_max precipitation_mm duree_ensoleillement annee
14/09/2013 4.18 13.6 21.5 49.5 154.2 2013
12/10/2013 4.19 11.7 17.8 43.4 96.5 2013
26/10/2013 4.22 11.7 17.8 43.4 96.5 2013
29/10/2013 4.30 11.7 17.8 43.4 96.5 2013
16/11/2013 4.39 5.7 10.0 69.7 43.7 2013
22/11/2013 4.48 5.7 10.0 69.7 43.7 2013

kable_styling permet un paramétrage de l’apparence du tableau, grâce au package kableExtra.

Ces données sont de type numériques, quantitatives, et continues.

Questions sciences

Quelles questions pourrait-on se poser avec des données comme ça ?

Est ce que le niveau de l’étiage a un lien avec les précipitations moyenne du mois, ou l’ensoleillement, ou les deux ? ou les temperatures ?

Y a t-il des interactions entre ces parametres ?
Est ce qu’il y a une saisonalité ? Est ce que les differences observees sont statistiquement significative ?
Quelles mesures supplementaires pourrait on faire pour consolider ses observations ?

Distribution des données

Dans un premier temps observons la distribution de ces valeurs. Suivant la distribution estimée, on pourra faire des tests paramétriques (distribution gaussienne de la population), ou des tests non-paramétriques (pas d’hypothèse quand à la distribution de la population)

L’histogramme des fréquences permet une visualisation de la distribution des données.

Histogramme

par(mfrow = c(3,2))

hist(as.vector(unlist(File1$etiage_m)), breaks=12,
     main = "Hauteur de l'etiage\nde la fontaine des capucins (2013-2015)",
     xlab = "hauteur (metres)", ylab = "Frequency", col = "aquamarine", prob=T)
curve(dnorm(x,mean(File1$etiage_m),sd(File1$etiage_m)), add=T,col="red")

hist(as.vector(unlist(File1$precipitation_mm)), breaks=10,
     main = "Precipitation moyenne mensuelle\na Paris (2013-2015)",
     xlab = "hauteur (mm)", ylab = "Frequency", col = "cadetblue1",prob=T)
curve(dnorm(x,mean(File1$precipitation_mm),sd(File1$precipitation_mm)), add=T,col="red")

hist(as.vector(unlist(File1$temp_mois_min)), breaks=12,
     main = "Temperature minimale \nmensuelle a Paris (2013-2015)",
     xlab = "temperature (celcius)", ylab = "Frequency", col = "chartreuse",prob=T)
curve(dnorm(x,mean(File1$temp_mois_min),sd(File1$temp_mois_min)), add=T,col="red")


hist(as.vector(unlist(File1$temp_mois_max)), breaks=12,
     main = "Temperature maximale \nmensuelle a Paris (2013-2015)",
     xlab = "temperature (celcius)", ylab = "Frequency", col = "chartreuse3",prob=T)
curve(dnorm(x,mean(File1$temp_mois_max),sd(File1$temp_mois_max)), add=T,col="red")


hist(as.vector(unlist(File1$duree_ensoleillement)), breaks=15,
     main = "duree d'ensoleillement moyen \nmensuelle a Paris (2013-2015)",
     xlab = "duree (minutes)", ylab = "Frequency", col = "chocolate2",prob=T)
curve(dnorm(x,mean(File1$duree_ensoleillement),sd(File1$duree_ensoleillement)), add=T,col="red")



par(mfrow = c(1,1))

par(mfrow = c(3,2)) permet de diviser la page en 2 lignes et 2 colonnes pour afficher les graphes.

Pour le moment cela ne donne pas des énormes informations, sauf que la distribution des données n’a pas l’air de suivre la loi normale (appelé aussi distribution gaussienne) qui donnerait un histogramme en forme de cloche avec la majorité des valeurs réuni au niveau de la moyenne). Il est important de savoir si les données suivent la loi normale, parceque l’on ne peut pas faire n’importe quels tests si ce n’est pas le cas.

test de Shapiro

Le test de Shapiro permet de dire si les variables suivent cette loi normale avec une p-value à comparer au seuil de probabilité d’erreur qu’on tolère (l’alpha).

# hauteur d'攼㸹tiage dans la fontaine (m)
shapiro.test(File1$etiage_m)

    Shapiro-Wilk normality test

data:  File1$etiage_m
W = 0.94281, p-value = 0.02729
# precipitation moyenne journaliere du mois (mm)
shapiro.test(File1$precipitation_mm)

    Shapiro-Wilk normality test

data:  File1$precipitation_mm
W = 0.91317, p-value = 0.002502
# temperature minimale du mois
shapiro.test(File1$temp_mois_min)

    Shapiro-Wilk normality test

data:  File1$temp_mois_min
W = 0.90317, p-value = 0.001191
# temperature maximale du mois
shapiro.test(File1$temp_mois_max)

    Shapiro-Wilk normality test

data:  File1$temp_mois_max
W = 0.90419, p-value = 0.001283
# duree d'ensoleillement
shapiro.test(File1$duree_ensoleillement)

    Shapiro-Wilk normality test

data:  File1$duree_ensoleillement
W = 0.93083, p-value = 0.01003

En test stat il y a une hypothèse nulle (h0) et une hypothèse alternative (h1).

h0 : les données suivent une loi normale
h1 : les données ne suivent pas une loi normale
alpha : par exemple 0.05 = accepter 5% de risque de se tromper

si la p-value (résultant du test de shapiro) est inférieure à un niveau alpha choisi (<0.05), alors l’hypothèse nulle est rejetée => il est improbable que ces données soient normalement distribuées => pas de distribution normale, donc on devra utiliser des tests non paramétriques.

si la p-value est supérieure au niveau alpha choisi (>0.05), alors on ne doit pas rejeter l’hypothèse nulle. On estime que la population a une distribution normale.

si on regarde les 5 resultats des tests shapiro : les p-value sont toutes <0.05 (donc on estime que cela ne suit pas la loi normale).

QQ plot

Le QQ plot permet également d’évaluer visuellement une distribution normale. Si les points sont alignés autour de la droite tracée la distribution de la variable étudiée est celle d’une loi normale.

par(mfrow = c(3,2))


qqnorm(File1$etiage_m, main ="hauteur d'etiage (m)")
qqline(File1$etiage_m,col="red")

qqnorm(File1$precipitation_mm, main ="precipitation (mm)")
qqline(File1$precipitation_mm,col="red")

qqnorm(File1$temp_mois_min, main ="temperature minimum (degre celcius)")
qqline(File1$temp_mois_min,col="red")

qqnorm(File1$temp_mois_max, main ="temperature maximum (degre celcius)")
qqline(File1$temp_mois_max,col="red")

qqnorm(File1$duree_ensoleillement, main ="duree d'ensoleillement (min)")
qqline(File1$duree_ensoleillement,col="red")


par(mfrow = c(1,1))

Même si les hypothèses de normalité sont rejetés pour l’échantillon, on peut choisir de prendre le risque de faire un test paramétrique, si on sait que la population est de distribution normale habituellement avec un plus grand nombre de données.

Correlations plot

voila, maintenant on veut voir les correlations entre nos paramètres. jour par jour.

Attention ! On confond souvent corrélation et regression :
La corrélation mesure l’intensité de la liaison entre des variables. La régression analyse la relation d’une variable par rapport à une ou plusieurs autres.

corr <- round(cor(File1[,2:6]), 1)

ggcorrplot(corr, method = "circle",  hc.order = TRUE, type = "lower")

On observe que les températures min et max (du meme jour) sont corrélées = les jours où les températures min sont basses les températures max sont basse aussi, et inversement.

La durée d’ensoleillement est corrélée avec les températures max = plus les jours sont long, plus les températures max sont fortes.

un peu corrélé : la durée d’ensoleillement avec les températures min = plus les jours sont longs plus les températures minimales sont élevé.

un peu corrélé : les précipitations avec avec les température min et max = quand il pleut beaucoup les températures min et max sont plus forte

pas corrélé : les précipitations avec la durée d’ensoleillement = A Paris il pleut toute l’année

un peu anticorrélé : le niveau d’étiage avec les précipitations en surface = plus il pleut, moins le niveau est haut (ah ! en voila un truc inatendu, va falloir réfléchir)

un peu anticorrélé : le niveau d’étiage avec la durée d’ensoleillement = plus les jours sont longs plus l’étiage est bas.

tres anticorrélé : le niveau d’étiage avec les températures min et max = plus le niveau d’étiage est haut plus les température sont faible.

La je pense qu’il faut se pencher sur la saisonnalité pour voir si les mécanismes ne sont pas influencé par ce parametre.

Serie temporelle

Nous allons faire un graphique de série temporelle.
D’abord il faut définir la colonne qui contient les dates.

Puis utiliser les fonctions graphiques de ggplot2 :
geom point : le nuage de points
geom_line : relier ces points
scale_x_date : permet de réduire les étiquettes des dates pour que ce soit un format lisible. stat_smooth : ajoute une ligne de tendance. ici j ai mis en paramètres la méthode “loess” qui est recommandée quand on a moins de 1000 observations.
ggtitle : le titre du graphe
xlab : pour le titre de l’axe des x
ylab : pour le titre de l’axe des y

Pour les couleurs, vous pouvez utiliser les noms des couleurs (ex: “red”, “ligthblue”, “chartreuse” etc…) ou par le code hexadécimal (ex: “#FFCC00”). Les codes se trouvent très facilement, taper “R color” dans un moteur de recherche.
Vous pouvez également utiliser des couleurs provenant du package RColorBrewer.

# d攼㸹finir la colone comme une variable de type date 
File1$date <- as.Date(File1$date, format = "%d/%m/%Y")

# v攼㸹rifier que c est bien consid攼㸹r攼㸹 comme des dates et que le format est correct
str(File1$date)
 Date[1:45], format: "2013-09-14" "2013-10-12" "2013-10-26" "2013-10-29" "2013-11-16" ...
# points reli攼㸹 : plot basique
ggplot(data = File1, aes(x = File1$date, y = File1$etiage_m)) + 
  geom_point() + 
  geom_line(color = "#00AFBB", size = 1, group = 1) +
  scale_x_date(date_labels = "%b/%Y")+ 
  stat_smooth(color = "#FC4E07", 
              fill = "#FC4E07",
              method = "loess")+
  ggtitle("Evolution de la hauteur d'etiage\nde la fontaine des Capucins")+
  xlab("Date")+
  ylab("hauteur en metres")

Ensuite j’aimerais afficher les saisons sur ce plot.
j’utilise la fonction geom_vline de ggplot2.
et geom_text pour ajouter le nom des saisons.

# points reli攼㸹 : plot basique
gp2 <- ggplot(data = File1, aes(x = File1$date, y = File1$etiage_m)) + 
  geom_point() + 
  geom_line(color = "#00AFBB", size = 1, group = 1) +
  scale_x_date(date_labels = "%b/%Y")+ 
  stat_smooth(color = "#FC4E07", 
              fill = "#FC4E07",
              method = "loess")+
  ggtitle("Evolution de la hauteur d'etiage\nde la fontaine des Capucins")+
  xlab("Date")+
  ylab("hauteur en metres")+
  geom_vline(xintercept = as.Date("2013-09-21"), linetype="dotted", color = "brown3", size=1)+
  geom_vline(xintercept = as.Date("2013-12-20"), linetype="dotted", color = "brown3", size=1)+
  geom_vline(xintercept = as.Date("2014-03-20"), linetype="dotted", color = "brown3", size=1)+
  geom_vline(xintercept = as.Date("2014-06-20"), linetype="dotted", color = "brown3", size=1)+
  geom_vline(xintercept = as.Date("2014-09-21"), linetype="dotted", color = "brown3", size=1)+
  geom_vline(xintercept = as.Date("2014-12-22"), linetype="dotted", color = "brown3", size=1)+

  geom_text(x=as.Date("2013-09-21")+45, y=4.6, label="Automne", color = "darkred")+
  geom_text(x=as.Date("2013-12-20")+45, y=4.6, label="Hivers", color = "darkred")+
  geom_text(x=as.Date("2014-03-20")+45, y=4.6, label="Printemps", color = "darkred")+
  geom_text(x=as.Date("2014-06-20")+45, y=4.6, label="Ete", color = "darkred")+
  geom_text(x=as.Date("2014-09-21")+45, y=4.6, label="Automne", color = "darkred")+
  geom_text(x=as.Date("2014-12-22")+45, y=4.6, label="Hivers", color = "darkred")  


# gp2

On aimerait avoir plus de données sur d’autres années pour vérifier ca mais il semblerait qu’il y ait une fluctuation cyclique de 50 cm du niveau de l’eau dans la fontaine d’étiage des Capucins, avec un pic en hivers et un creux en été.

Maintenant je voudrais savoir si cette fluctuation entre saison est statistiquement significative. Pour ca je vais devoir utiliser un test statistique. Comme on a vu au début que les données de hauteur d’étiage n’avaient pas une distribution suivant la loi normale (histogrammes, et shapiro test), on va effectuer un test non paramétrique appelé Kruskal Wallis. Mais d’abord on va devoir ajouter une colonne “saison” à notre tableau. D’ailleurs, il aurait d’ailleurs été plus malin de faire cela avant pour éviter de taper manuellement les positions des lignes et texte dans le graphique précédent. Mais comme ca vous savez bidouiller aussi.

Meme chose pour les temperatures. comme les deux températures sont des colonnes differentes je remet l’aes avec les définitions des x et y pour la deuxieme, et elle s’ajoutte au meme plot en la mettant dans un objet gpl auquel j’additione la suite.

# points reli攼㸹 : plot basique
gpl <- ggplot(data = File1, aes(x = File1$date, y = File1$temp_mois_min)) + 
  geom_point(color="aquamarine4") + 
  geom_line(color = "aquamarine3", size = 1, group = 1) +
  scale_x_date(date_labels = "%b/%Y")+ 
  stat_smooth(color = "chartreuse4", 
              fill = "chartreuse",
              method = "loess")
  
gpl <- gpl+ geom_point(data =File1, aes(x = File1$date, y = File1$temp_mois_max)) + 
  geom_point(color="aquamarine4") +
  geom_line(data =File1, aes(x = File1$date, y = File1$temp_mois_max), 
            color = "coral4", 
            size = 1, group = 1) +
  stat_smooth(data =File1, aes(x = File1$date, y = File1$temp_mois_max), color = "#FC4E07", 
              fill = "coral",
              method = "loess")

gpl <- gpl +  ggtitle("Temperature min et max en surface")+
  xlab("Date")+
  ylab("degre celcius")+
  geom_vline(xintercept = as.Date("2013-09-21"), linetype="dotted", color = "brown3", size=1)+
  geom_vline(xintercept = as.Date("2013-12-20"), linetype="dotted", color = "brown3", size=1)+
  geom_vline(xintercept = as.Date("2014-03-20"), linetype="dotted", color = "brown3", size=1)+
  geom_vline(xintercept = as.Date("2014-06-20"), linetype="dotted", color = "brown3", size=1)+
  geom_vline(xintercept = as.Date("2014-09-21"), linetype="dotted", color = "brown3", size=1)+
  geom_vline(xintercept = as.Date("2014-12-22"), linetype="dotted", color = "brown3", size=1)+

  geom_text(x=as.Date("2013-09-21")+45, y=0, label="Automne", color = "darkred")+
  geom_text(x=as.Date("2013-12-20")+45, y=0, label="Hivers", color = "darkred")+
  geom_text(x=as.Date("2014-03-20")+45, y=0, label="Printemps", color = "darkred")+
  geom_text(x=as.Date("2014-06-20")+45, y=0, label="Ete", color = "darkred")+
  geom_text(x=as.Date("2014-09-21")+45, y=0, label="Automne", color = "darkred")+
  geom_text(x=as.Date("2014-12-22")+45, y=0, label="Hivers", color = "darkred")  


# gpl

Meme chose pour les precipitations en surface.

# points reli攼㸹 : plot basique
gp1 <- ggplot(data = File1, aes(x = File1$date, y = File1$precipitation_mm)) + 
  geom_point(color="blue") + 
  geom_line(color = "blue4", size = 1, group = 1) +
  scale_x_date(date_labels = "%b/%Y")+ 
  stat_smooth(color = "darkblue", 
              fill = "cornflowerblue",
              method = "loess")+
  ggtitle("Precipitations en surface")+
  xlab("Date")+
  ylab("precipitation en mm")+
  geom_vline(xintercept = as.Date("2013-09-21"), linetype="dotted", color = "brown3", size=1)+
  geom_vline(xintercept = as.Date("2013-12-20"), linetype="dotted", color = "brown3", size=1)+
  geom_vline(xintercept = as.Date("2014-03-20"), linetype="dotted", color = "brown3", size=1)+
  geom_vline(xintercept = as.Date("2014-06-20"), linetype="dotted", color = "brown3", size=1)+
  geom_vline(xintercept = as.Date("2014-09-21"), linetype="dotted", color = "brown3", size=1)+
  geom_vline(xintercept = as.Date("2014-12-22"), linetype="dotted", color = "brown3", size=1)+

  geom_text(x=as.Date("2013-09-21")+45, y=4.6, label="Automne", color = "darkred")+
  geom_text(x=as.Date("2013-12-20")+45, y=4.6, label="Hivers", color = "darkred")+
  geom_text(x=as.Date("2014-03-20")+45, y=4.6, label="Printemps", color = "darkred")+
  geom_text(x=as.Date("2014-06-20")+45, y=4.6, label="Ete", color = "darkred")+
  geom_text(x=as.Date("2014-09-21")+45, y=4.6, label="Automne", color = "darkred")+
  geom_text(x=as.Date("2014-12-22")+45, y=4.6, label="Hivers", color = "darkred")  

# gp1

juste pour le plaisir de les afficher a la suite

gpl

gp1

gp2

Ajouter une colonne “saison”

Bon pour ajouter la colone saison y’a forcément plus simple, mais ca reste une solution qui marche.

# je separe le tableau par saison avec subset
automne13 <- subset(File1, date > as.Date("2013-09-21") & date < as.Date("2013-12-22"))
hivers13 <- subset(File1, date > as.Date("2013-12-22") & date < as.Date("2014-03-20"))
printemps14 <- subset(File1, date > as.Date("2014-03-20") & date < as.Date("2014-06-20"))
ete14  <- subset(File1, date > as.Date("2014-06-20") & date < as.Date("2014-09-21"))
automne14  <- subset(File1, date > as.Date("2014-09-21") & date < as.Date("2014-12-22"))
hivers14 <- subset(File1, date > as.Date("2014-12-22") & date < as.Date("2015-03-20"))

# rbind permet de rassemble 2 tableaux ayant les memes colonnes, ici pour rassembler les saisons d annees differentes.
automne <- rbind(automne13, automne14)
hivers <- rbind(hivers13, hivers14)
printemps <- printemps14
ete <- ete14

# ajouter la colonne saison dans chaque sous tableau
automne$saison <- "automne"
hivers$saison <- "hivers"
printemps$saison <- "printemps"
ete$saison <- "ete"


# reassembler le tableau final

File0 <- rbind(automne, hivers, printemps, ete)
# View(File0)

Analyse en composante principale

** analyse factorielle **

Cette analyse permet de visualiser des données décrites avec plusieurs variables quantitatives, et de voir comment ces échantillons (= les jours de mesures de hauteurs d’eau) sont influencés par ces variables (la temperature min, max, précipitation, duree d’ensoleillement, et hauteur d’etiage),

res.pca <- prcomp(File0[, 2:6],  scale = TRUE)


# ACP indiv
p <- fviz_pca_ind(res.pca, label="none", habillage=File0$saison,
             addEllipses=TRUE, ellipse.level=0.95)
print(p)

#Graphique des variables
fviz_pca_var(res.pca, col.var="steelblue")+
 theme_minimal()

# indiv et variables
fviz_pca_biplot(res.pca, label="var", habillage=File0$saison,
               addEllipses=TRUE, ellipse.level=0.95)

On observe que 71,% de la variation est expliqué par l’axe 1 (Dim1) et que sur cet axes la variable de hauteur d’étiage est anti-corrélée avec toutes les autres. Les groupes de saisons (les élipses colorées) ne sont pas très distinct.

Tests de comparaison

Ici voici des exemples de tests proposant de comparer des groupes d’échantillons (variable qualitatives ou facteur, ici, la saison), et d’estimer leur effet sur une variable quantitative (ex. : la hauteur d’étiage…)

Test Kruskal Wallis

** test non-paramétrique ** A faire lorsque les données ne sont pas distribuées selon la loi normale.

Ce test répond à la question y a t-il une difference significative des hauteurs d’étiage de la fontaine des Capucins entre les differentes saisons ?

Oui : si la p-value du test est < 0.05 (à un seuil de 5% d’erreur acceptée)

# test hauteur 攼㸹tiage
testKW <- kruskal.test(etiage_m~saison , data = File0)
print(testKW)

    Kruskal-Wallis rank sum test

data:  etiage_m by saison
Kruskal-Wallis chi-squared = 22.727, df = 3, p-value = 4.604e-05
# si besoin de la pvalue
# pval <- as.character(testKW[3]$p.value)

La p-value < 0.05 , donc il y a bien une difference significative entre les saisons pour la hauteur d’etiage (avec 5% d’erreur acceptée).

Test ANOVA

** test paramétrique ** * pour distribution de la population assumée normale*

ANOVA : analyse de la variance

Avant de faire ce test il faut vérifier que la variance entre les groupes est homogène. si p-value du test bartlett > 0.05 : les variances des groupes sont homogenes et on peut faire le test ANOVA.

# tester l homoscedasticite
bartlett.test(etiage_m ~ saison, data = File0)

    Bartlett test of homogeneity of variances

data:  etiage_m by saison
Bartlett's K-squared = 16.905, df = 3, p-value = 0.0007391
# test ANOVA
testAOV <- aov(etiage_m~saison , data = File0)
AOVi <- summary(testAOV)
print(AOVi)
            Df Sum Sq Mean Sq F value   Pr(>F)    
saison       3 0.5561 0.18536   14.62 1.41e-06 ***
Residuals   40 0.5070 0.01267                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# si besoin de la pvalue
# pval2 <- as.character(AOVi[[1]]$Pr[1])
# print(pval2)

test bartlett

la p-value < 0.05 donc on ne devrait pas faire une ANOVA, mais le test Kruskal Wallis etait plus approprié.
Encore une fois ce probleme serait probablement corrigé si on avait plus de mesures.

test ANOVA (juste pour la démo, je l’ai quand même fait)
La p-value < 0.05 , donc il y a bien une difference significative entre les saisons pour la hauteur d’etiage (avec 5% d’erreur acceptée).

Post-hoc TukeyHSD

Maintenant on veut savoir quel sont les saisons qui sont differentes les unes des autres ? pour ca il faut faire un test post-hoc. Ici je propose le test TukeyHSD.

tuk <- TukeyHSD(testAOV)
kable(tuk$saison) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

diff lwr upr p adj
ete-automne -0.1208788 -0.2406674 -0.0010902 0.0472330
hivers-automne 0.1850303 0.0652417 0.3048189 0.0009641
printemps-automne -0.0549048 -0.1930346 0.0832251 0.7122280
hivers-ete 0.3059091 0.1772355 0.4345827 0.0000008
printemps-ete 0.0659740 -0.0799281 0.2118762 0.6229745
printemps-hivers -0.2399351 -0.3858372 -0.0940329 0.0004291

Boxplot

Les boxplot ou boites à moustaches permettent de visualiser la dispersion des valeurs pour chaque classe (les saisons) et leur médiane (le trait noir central qui sépare 50% des valeurs).

bplot <- ggplot(File0,
  aes(x=saison,
      y=etiage_m,
      fill=(factor(saison))))

bplot <- bplot + geom_boxplot(width=.8,
  outlier.size=1.25,
  outlier.shape=21)+ ggtitle("Saisons et hauteur d'etiage")+ 
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

bplot <- bplot + geom_point(position=position_jitterdodge(dodge.width=0.7), size=1)

bplot

On constate que les hauteurs d’étiage en automne sont plus aléatoire que celle des 3 autres saisons.

On aimerait pouvoir voir les binomes de saisons significativement differentes pour la hauteur d’étiage identifiée par le test TukeyHSD.

Cette fois je vais utiliser le visuel violin plot

Violin plot avec p-val

Dans la liste créée my_comparison, je mets que les binomes qui sont significatif sur le tableau de résultat du test TukeyHSD (p-value < 0.05). les p-values significatives sont affichées au dessus des violin plot.

##boxplot summary genotype
my_comparisons <- list(c("automne", "hivers"), c("ete", "hivers"), c("printemps", "hivers"))

p <- ggplot(File0, aes(x = saison, y = etiage_m, fill = factor(saison)))

p <- p+ geom_violin(width=.8)+ 
  ggtitle("Saisonnalite de la hauteur d'etiage") + 
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  geom_dotplot(binaxis='y', stackdir='center', dotsize=0.3)

p <- p + geom_boxplot(width=0.1, fill="white")
p + stat_compare_means( comparisons = my_comparisons, test = "kruskal.test", 
                        label.y = c(150) +
theme(legend.position = "none"))

Tricoter son Rmd

Une fois que vous avez fini votre script Rmd, vérifié que tout les chunk fonctionnent un par un, vous pouvez l’exécuter en entier pour produire le fichier final (Html dans mon cas). Il suffit de cliquer sur Knit / Knit to html

ici :

chunkadd

chunkadd

Exporter le tableau

Comme le tableau a été un peu modifié (ajout d’une colonne) on pourrait avoir envie de le sauver sous sa nouvelle forme en .csv pour le réutiliser plus tard dans un autre script… par exemple

write.table(File0,
            file=file.path(dir.base,
                "K_tab_v02.csv"),
            sep=';', 
            col.names = TRUE, 
            row.names = FALSE, 
            quote = FALSE)

Conclusion

On vous pose souvent la question est ce que le niveau d’eau varie beaucoup dans les catacombes quand il pleut ?
Vous pourrez répondre que il varie oui, significativement en hivers par rapport au autres saisons, d’environ un bon 20cm plus haut que les autres saisons. Non ce n’est pas directement corrélés aux précipitations, mais plutot au température plus basses. Ou devrait-on dire que le niveau est maximale en hivers et baisse significativement lors des autres saisons.

Hypothèse d’explication : Les températures plus importantes et épisodes de sécheresse doivent engendrer une forte évaporation des apports potentiels en eaux depuis la surface avant que celles-ci ne s’infiltre pour réalimenter le sous-sol. Donc en réponse version courte : Non le niveau ne monte ni beaucoup ni brusquement en cas de forte pluie, mais par contre le niveau baisse quand il fait chaud. Et ce phénomène est saisonnier, et progressif.

Mais attention ces conclusions ne fonctionnent que sur cet échantillon de données que l’on possède et qui est très restreint, une année et quelque, mesures une fois toute les semaines environ. Il faudrait trouver de nouvelles données et vérifier ce qui se passe. Plus de jour de mesures, à des fréquences plus courte, d’autres années ayant des météos différentes, différentes stations de mesure.

Est ce qu’ils n’ont pas été un peu pessimistes de faire des échelles d’étiage si hautes ?

Qui est motivés pour prendre des mesures régulièrement au bar des rats, à Cubes et aux Chartreux pour vérifier si ca se confirme ?