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 ?

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…

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.

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 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"

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 !

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 version basique
# kable(head(File1))

# la version parametrable
kable(head(File1)) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

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.

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")

# 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.

# 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

# 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

gpl

gp1

gp2

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.

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 ?