Dans un précédent article j’abordais, d’un point de vue théorique, les méthodes qui peuvent être employées pour évaluer la liaison entre deux variables numériques continues. Je disais qu’il en existe deux principales : la régression linéaire simple qui permet d’évaluer s’il existe un lien linéaire entre les deux variables, et la corrélation de Spearman qui permet d’évaluer l’existence d’un lien moins spécifique : la monotonie.

Ces deux approches sont décrites en détail dans [ce premier article. Vous y trouverez également leurs conditions d’applications, ainsi que leur points de similitude et de divergence.Dans ce second volet qui est un tutoriel R, je vais vous montrer comment utiliser la régression linéaire simple et la corrélation de Spearman avec le logiciel R.

Pour illustrer cet article j’ai choisi d’utiliser un jeu de données issu d’une expérimentation de terrain qu’un doctorant avec lequel je travaille a du analyser. Les données de terrains ne sont généralement pas des cas idéaux, en partie parce que les sources de variabilité sont moins bien contrôlées qu’en laboratoire. De ce fait, elles reflètent davantage les difficultés auxquelles on peut être confronté lorsqu’on souhaite évaluer si deux variables numériques continues sont liées.

1. Les données

Le domaine dans lequel ces données ont été acquises est très spécifique. Alors, dans un soucis principalement de compréhension, nous allons imaginer qu’une des variables est le pourcentage de surfaces cultivées au niveau d’un site de prélèvement, et que la seconde est la concentration d’un élément chimique quelconque mesurée sur ce même site. Au total 28 sites ont été étudiés.

mydata <-structure(list(conc = c(0.72, 0.72, 0.88, 0.99, 1.2, 1.303538175,
0.93, 0.97, 1.54, 1.47, 1.61, 2.5, 1.29, 1.17, 1.06, 1.9, 1.88,
2.06, 1.76, 2.37, 2.13, 2.12, 2.11, 1.81, 7.784433, 5.1, 1.95,
3.49), pSurface_cult = c(6.269481783, 0, 0, 100, 77.55458855,
0, 22.04150899, 40.40614454, 40.60597188, 22.04150899, 22.04150899,
77.55458855, 40.40614454, 100, 52.58843992, 34.12077483, 98.94265481,
48.12059561, 67.60579421, 22.04150899, 5.024702025, 67.60579421,
42.14948675, 52.58843992, 60.99253217, 48.12059561, 100, 100)), class = "data.frame", row.names = c(NA,
-28L), .Names = c("conc", "pSurface_cult"))

2. Premières visualisations de la liaison entre deux variables

Lorsqu’on souhaite étudier la liaison entre deux variables numériques continues, la première chose à faire est de visualiser les données à l’aide d’un scatter-plot. Il s’agit tout d’abord d’évaluer la forme de la relation entre les deux variables.

En effet, pour rappel, la régression linéaire simple, ne peut être employée que si la relation entre les deux variables est globalement de forme linéaire. De même, la corrélation de Spearman ne peut être employée que si la relation est globalement monotone. L’appréciation de la forme de la relation (linéaire ou monotone) se fait uniquement de façon visuelle.

Le scatter plot permet également d’évaluer le sens de la relation : est ce que les deux variables varient dans le même sens ou en en opposé?

Il permet enfin de mettre en évidence la présence éventuelle d’outliers. Leur détection est importante car ils peuvent entraîner un défaut de normalité des résidus, et de la statistique T employée pour tester la significativité de la pente de la régression. Et au final, conduire à une conclusion erronée au sujet de la relation entre les deux variables.

Je réalise quasiment toutes mes représentations graphiques avec le package ggplot2. Ce package nécessite un effort d’apprentissage mais ensuite il est extrêmement pratique et performant. A tel point qu’il est devenu incontournable dans le domaine de l’analyse de données. Pour ceux qui ne connaissent pas encore ce package, vous trouverez un article d’introduction à la visualisation sous R avec le package ggplot2, en cliquant ici, et sa feuille d’aide (cheat sheet) là.Dans un premier scatterplot, on peut ajouter une courbe de régression locale (Loess)  pour aider à visualiser la forme globale de la relation entre les deux variables.

library(ggplot2)

ggplot(mydata, aes(y=conc, x=pSurface_cult))+
geom_point()+
geom_smooth(colour="blue", method="loess", se=FALSE)

Ce plot montre une relation plutôt biphasique entre les deux variables avec une première phase croissante , suivie d’une phase de décroissance. Néanmoins, la deuxième phase semble fortement influencée par deux points extrêmes. La monotonie semble donc contrariée principalement par 2 points sur 28, soit environ 7% des données. De ce fait, on peut penser que l’hypothèse d’une relation globalement monotone entre les deux variables n’est pas complètement absurde.

Dans un deuxième temps, on peut ajouter une droite de régression linéaire sur le scatter plot, pour évaluer si l’hypothèse d’une tendance linéaire entre les deux variables est acceptable.

ggplot(mydata, aes(y=conc, x=pSurface_cult))+
geom_point()+
geom_smooth(colour="red", method="lm", se=FALSE)

Ici encore, l’hypothèse d’une relation globalement linéaire entre les deux variables ne semble pas absurde. Néanmoins, l’influence potentielle des deux points extrêmes sur le test de la significativité de le pente (autrement dit sur le test de l’égalité à 0 de la pente) devra être explorée.

3. Evaluation d’une liaison linéaire par régression linéaire simple

Comme expliqué en détail dans le premier article, la régression linéaire permet de quantifier et d’évaluer la significativité d’une liaison linéaire entre deux variables numériques continues.Nous allons considérer ici que :

•  la concentration (“conc”) est la variable réponse
• le pourcentage de surfaces cultivées (“pSurface_cult”) est la variable prédictive.

3.1 Ajustement du modèle linéaire simple

L’ajustement du modèle linéaire par la méthode des moindres carrés ordinaires est réalisé par la fonction lm.

mod_lm <-lm(conc ~ pSurface_cult,data=mydata)

3.2 Evaluation des hypothèses

Le test de significativité de la pente est valide si :

• les réponses sont indépendantes,
• les résidus de la régression sont distribués selon une loi Normale
• les résidus de la régression sont homogènes

L’hypothèse de l’indépendance des réponses est, à priori, validée par le protocole, nous allons considéré que c’est le cas ici. Les deux autres hypothèses en revanche doivent être évaluées.

3.2.1 Evaluation de l’hypothèse de normalité des résidus

L’hypothèse de normalité des résidus est explorée visuellement et à l’aide d’un test statistique.Plusieurs méthodes visuelles et plusieurs test peuvent être employés mais les plus courantes sont le QQplot et le test deShapiro_Wilk.

plot(mod_lm, 2)

Pour accepter la normalité des résidus, ceux si doivent être alignés le long de la droite tracée en pointillés. C’est globalement le cas ici, à l’exception des deux points extrêmes 25 et 26 déjà mis en évidence sur les scatter plots.

L’hypothèse de normalité des résidus est rejetée lorsque le QQ plot a une forme de banane comme ci dessous :

shapiro.test(residuals(mod_lm))
##
## Shapiro-Wilk normality test
##
## data: residuals(mod_lm)
## W = 0.68821, p-value = 1.95e-06

Pour accepter la normalité, la p-value du test de Shapiro Wilk doit être supérieure au seuil de 0.05. Dans notre cas, la p-value est largement inférieure au seuil de significativité, l’hypothèse est donc rejetée. Ce rejet est vraisemblablement du aux deux points extrêmes.

 

3.2.2 Evaluation de l’hypothèse d’homogénéité des résidus

L’hypothèse d’homogénéité des résidus peut également être évaluée à la fois visuellement et par l’utilisation d’un test statistique.La méthode visuelle consiste à réaliser un plot de la racine carrée des résidus standardisés vs. les valeurs ajustées par le modèle. Cela peut facilement être réalisé avec la commande suivante :

plot(mod_lm, 3)

L’hypothèse d’homogénéité est acceptée si les points sont globalement distribués dans un rectangle ; c’est le cas ici.

Elle est en revanche rejetée si les points sont distribués selon un motif systématique comme en entonnoir par exemple (c’est à dire un triangle avec un sommet à droite et une base à gauche, ou l’inverse).

Le test statistique classiquement employé pour évaluer l’homogénéité des résidus est un test de score, parfois appelé test de Breusch-Pagan. Il est disponible dans le package `car` avec la fonction `ncvTest`, ou encore sous une version “studentisée” (plus robuste selon certains), dans le package `lmtest` avec la fonction `bptest`.

library(car)
ncvTest(mod_lm)

## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 1.639657 Df = 1 p = 0.2003726

library(lmtest)
bptest(mod_lm)

##
## studentized Breusch-Pagan test
##
## data: mod_lm
## BP = 0.32891, df = 1, p-value = 0.5663

Le plot ne montre pas de défaut manifeste d’homogénéité, et les p-values des deux versions du test de Breush Pagan est supérieures à 0.05. Au final l’hypothèse d’homogénéité des résidus est acceptée.

L’hypothèse de normalité ayant été rejetée par le test de Shapiro Wilk, une transformation log10 peut être appliquée à la variable réponse, (ici “conc”) pour essayer de l’améliorer.

3.3 Utilisation d’une transformation log10 des réponses

3.3.1 Ajustement du modèle

mod_lm2 <-lm(log10(conc) ~ pSurface_cult, data=mydata)

 

3.3.2 Evaluation des hypothèses de normalité et d’homogénéité des résidus

plot(mod_lm2,2)

 

shapiro.test(residuals(mod_lm2))
## Shapiro-Wilk normality test
##
## data: residuals(mod_lm2)
## W = 0.92771, p-value = 0.05399

Sur le QQ plot les points 25 et 26 semblent moins extrêmes et la pvalue du test de Shapiro-Wilk est cette fois supérieure (mais tout juste) à 0.05. Au final, les résidus souffrent d’un défaut de normalité, mais celui ci reste dans la limite de l’acceptable. Il faudra néanmoins rester vigilant sur la validité des résultats.

plot(mod_lm2,3)

ncvTest(mod_lm2)

## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.3792163 Df = 1 p = 0.538023

bptest(mod_lm2)

##
## studentized Breusch-Pagan test
##
## data: mod_lm2
## BP = 0.24846, df = 1, p-value = 0.6182

Le plot et les tests ne mettent pas en évidence de défaut d’homogénéité des résidus. L’hypothèse d’homogénéité des résidus est donc acceptée.

3.4 Identification des données “pas comme les autres”

Il est toujours intéressant de détecter si certaines données se comportent de manière différente des autres. C’est à dire si certaines données peuvent être considérées comme des outliers et/ ou comme ayant une influence importante sur les paramètres de la régression (la pente et l’ordonnée à l’origine).Le package car dispose d’une fonction très utile pour détecter ces données différentes, il s’agit de la fonction influenceIndexPlot. Cette fonction renvoie 4 graphs :

• le premier (en partant du bas), celui des hat value, reflète l’effet de levier (ou poids) de chaque donnée sur sa propre estimation. Une donnée est considérée comme atypique lorsque cette valeur est <0.5.
• le second plot celui des p-value de Bonferroni permet de mettre en évidence les outliers. Est considérée comme outlier une donnée ayant une p-value inférieure à 0.05.
• le troisième plot, celui des résidus studentizés permet également de mettre en évidence les outliers
• le dernier plot, celui des distance de Cook permet d’évaluer l’influence des données sur les paramètres de régression. La distance de Cook mesure le changement dans l’estimation des paramètres de régression lorsque la donnée n’est pas prise en compte par les moindres carrés. Plus la distance est élevée, plus la modification des paramètres de régression est importante. Le seuil de 1 est couramment utilisé pour considéré qu’une donnée est très influente.

Vous trouverez plus de détail sur l’effet de levier et les distances de Cook là.

influenceIndexPlot(mod_lm2)

La donnée 25 est mise en évidence par 3 des 4 plots. Les deuxième et troisième graph (en partant du bas) la désigne comme potentiel outlier. De son côté, le plot des distances de Cook montre que son influence est parmi les deux plus fortes. Néanmoins elle est inférieure à 1, ce qui nous amène à penser que son influence sur les paramètres du modèle n’est pas vraiment problématique. La pvalue de Bonferroni (plot 2) peut être obtenue avec la fonction `outlierTest(nom_du_modele)`.

outlierTest(mod_lm2)
##
## No Studentized residuals with Bonferonni p< 0.05 
## Largest |rstudent|: 
## rstudent unadjusted p-value Bonferonni p 
## 25 3.398425 0.0022757 0.06372 

La pvalue ajustée par la méthode de Bonferonni est très proche du seuil de 0.05. A mon sens on peut considérer la donnée 25 comme une outlier. Il sera donc intéressant de prendre connaissance des résultats de la régression linéaire avec et sans cette donnée. En revanche, la donnée 26 ne semble pas poser de problèmes particuliers.

3.5 Résultats

Les hypothèses de normalité et d’homogénéité étant acceptées, les résultats de la régression sont valides, nous pouvons donc les consulter à l’aide de la fonction summary(nom_du_modele)

summary(mod_lm2) ## 
## Call: 
## lm(formula = log10(conc) ~ pSurface_cult, data = mydata) 
## 
## Residuals: ## Min 1Q Median 3Q Max 
## -0.34015 -0.19422 -0.00441 0.09968 0.64357 
## 
## Coefficients: 
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.109859 0.077345 1.420 0.167
## pSurface_cult 0.002259 0.001330 1.699 0.101
##
## Residual standard error: 0.2293 on 26 degrees of freedom
## Multiple R-squared: 0.09991, Adjusted R-squared: 0.06529
## F-statistic: 2.886 on 1 and 26 DF, p-value: 0.1013

La partie `Residuals` des résultats permet d’évaluer rapidement la normalité des résidus. Lorsque les résidus sont distribués selon un loi Normale, la médiane doit être autour de 0 (c’est le cas ici), et les valeurs absolues de Q1 (premier quartile) et Q3 (troisième quartile) doivent être proches. Ce n’est pas tout à fait le cas ici, mais nous savions déjà que les résidus souffrent d’un défaut de normalité.

La première ligne de la partie coefficients concerne l’ordonnée à l’origine, alors que la seconde ligne concerne la pente. Concernant les colonnes :

• la première colonne rapporte l’estimation des coefficients des paramètres,
• la seconde colonne l’estimation de leur erreur standard,
• la troisième colonne est la statistique T
• la dernière colonne rapporte la p-value du test évaluant l’égalité à 0 des coefficients.

La p-value du test de significativité de la pente est égale à 0.101. De ce fait, si on applique strictement le seuil de significativité classiquement fixé à 0.05, alors l’hypothèse de l’égalité à 0 de la pente est acceptée. Autrement dit, aucun lien linéaire significatif entre les deux variable n’est mis en évidence.

Si on souhaite caractériser l’évolution de la concentration en fonction du pourcentage de surface cultivée, les résultats nous montrent qu’en moyenne le log10 de la concentration augmente d’environ 2.26.10^-3 lorsque la surface cultivée augmente d’un pour cent.

Regardons à présent ce que deviennent ces résultats lorsque la donnée 25 est retirée. Pour cela, on refait tourner la régression sans la donnée, et on utilise la fonction `compareCoefs` du package `car`.

mod_lm3 <-lm(log10(conc) ~ pSurface_cult, data=mydata[-25,]) 
compareCoefs(mod_lm2, mod_lm3) 
## 
## Call: 
## 1: lm(formula = log10(conc) ~ pSurface_cult, data = mydata) 
## 2: lm(formula = log10(conc) ~ pSurface_cult, data = mydata[-25, ]) 
## Est. 1 SE 1 Est. 2 SE 2 
## (Intercept) 0.10986 0.07734 0.09983 0.06530 
## pSurface_cult 0.00226 0.00133 0.00197 0.00112

Les résultats nous montrent que les estimations des paramètres et de leur erreurs standards de ne sont effectivement pas très différents avec et sans la donnée 25. Qu’en est il de la p-value du test de l’égalité à 0 de la pente ?

summary(mod_lm3) 
## 
## Call: 
## lm(formula = log10(conc) ~ pSurface_cult, data = mydata[-25, ]) 
## 
## Residuals: 
## Min 1Q Median 3Q Max 
## -0.30117 -0.17409 0.01251 0.11548 0.51295 
## 
## Coefficients: ## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.099830 0.065301 1.529 0.1389
## pSurface_cult 0.001970 0.001125 1.751 0.0922 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1934 on 25 degrees of freedom
## Multiple R-squared: 0.1093, Adjusted R-squared: 0.07362
## F-statistic: 3.066 on 1 and 25 DF, p-value: 0.0922

La pvalue du test est à présent de 0.09.

Lorsque la pvalue d’un test est comprise dans l’intervalle [0.05 ;0.1] certains considèrent qu’elle est presque significative. Et de ce fait, qu’elle constitue une faible preuve (weak evidence ) de la relation linéaire. Dans cette situation le terme de “tendance” est généralement employé, comme par exemple : “les résultats mettent en évidence une tendance linéaire entre les deux variables”

D’autres considèrent que les pvalues comprises dans l’intervalle [0.03 ; 0.07] sont dans une zone “grise” qui ne permet pas de conclure sur l’acceptation ou le rejet de l’hypothèse H0 (ici l’égalité à 0 de la pente).

La mise en évidence d’un outlier doit inciter à vérifier cette donnée. Ici, nous n’avons pas de raison de penser qu’il s’agit d’une erreur, il n’y a donc pas de raison de prendre en compte les résultats de cette régression plutôt que ceux du modèle “mod_lm2”.

3.6 Discussion

Au final, il n’est pas facile d’avoir un avis tranché sur l’existence ou non d’une relation linéaire significative entre ces deux variables compte tenu :

• qu’une transformation log de la concentration à du être employée pour satisfaire l’ hypothèse de normalité des résidus. De ce fait, la relation qui est étudiée n’est plus celle entre la concentration et le pourcentage de surfaces cultivées mais entre le log10 de la concentration et le pourcentage de surfaces cultivées
• que malgré la transformation log de la réponse, la normalité des résidus reste médiocre
• que la taille de l’échantillon est assez faible : 28
• qu’une donnée est outlier
• et que la pvalue est égale à 0.1.

Si seule l’évaluation du lien linéaire avait eu un intérêt ici, il aurait été difficile de conclure sereinement. Dans le cas de ces données, la question posée était seulement d’identifier, si lorsque l’une des variables augmente, la seconde augmente aussi, ou au contraire diminue. L’évaluation de la monotonie est donc une autre approche possible.

3.7 Représentation finale de la régression

Lorsque la pente est significative, on réalise généralement une représentation graphique pour illustrer la relation linéaire qui lie les deux variables. Cette visualisation doit être en adéquation avec la régression : si une transformation a été utilisée dans le modèle de régression (ici le log10 de la réponse), elle doit également figurer dans la représentation graphique.

En plus des données observées et de la droite de régression, il est également de coutume de représenter l’intervalle de confiance à 95% de la pente. La fonction `geom_smooth` du package `ggplot2` permet de le faire facilement car c’est l’option par défaut. On peut également ajouter l’équation de la droite, à l’aide de la fonction `annotate`.

ggplot(mydata, aes(y=log10(conc), x=pSurface_cult))+
geom_point()+
geom_smooth(colour="red", method="lm") +
ylab("Concentration ")+
annotate("text", x = 25, y = 0.8, label = "log10(conc) = 0.11 + 0.00226 * pSurf_ cult\n (pval=0.1)")

 

4. Evaluation du lien monotone par la corrélation de Spearman

Nous avons vu précédemment, grâce aux premiers scatter plots, que l’hypothèse d’une relation globalement monotone entre les deux variables n’était pas aberrante. Nous allons donc caractériser la force de cette relation et évaluer sa significativité en utilisant la corrélation de Spearman.

 

4.1 Calcul du coefficient de corrélation de Spearman et évaluation de sa signficativité

La fonction cor.test du package stats (chargé par défaut à chaque ouverture de session R) permet de calculer le coefficient de corrélation de Spearman et de tester sa significativité.Cette fonction dispose de plusieurs arguments qui permettent, entre autre, de :

• sélectionner la méthode de corrélation souhaitée (Pearson Spearman ou Kendall)
• sélectionner le type d’hypothèse alternative souhaitée : uni ou bilatérale.
• préciser le type de p-value souhaitée : exacte ou par approximation normale.

Pour plus de détail, consulter l’aide sur cette fonction.

Comme nous n’avons pas d’à priori sur le sens de la relation, nous utilisons une hypothèse alternative bilatérale (option par défaut).

cor.test(mydata$conc, mydata$pSurface_cult, method="spearman")

## Warning in cor.test.default(mydata$conc, mydata$pSurface_cult, method =
## "spearman"): Cannot compute exact p-value with ties

##
## Spearman's rank correlation rho
##
## data: mydata$conc and mydata$pSurface_cult
## S = 2404.9, p-value = 0.07499
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.3418547

Le coefficient de corrélation est de 0.34. Cela peut paraître faible mais il n’existe pas de seuil universel (c’est à dire applicable à tous les domaines d’études)à partir duquel on pourrait considérer que larelation est forte. La p-value est égale à 0.0749, elle est proche de la significativité.

Le warning renvoyé par le logiciel R nous indique seulement que les données comportent des ex-aequos (au niveau des pourcentages de surfaces cultivées), et que cela empêche le calcule de la p-value selon une méthode exacte. La p-value est donc dérivée d’une approximation normale.

Compte tenu que la pvalue se situe dans l’intervalle [0.05, 0.1], il peut être intéressant de d’utiliser, de façon complémentaire, le coefficcient de corrélation de Kendall.

cor.test(mydata$conc, mydata$pSurface_cult, method="kendall")

## Warning in cor.test.default(mydata$conc, mydata$pSurface_cult, method =
## "kendall"): Cannot compute exact p-value with ties
##
## Kendall's rank correlation tau
##
## data: mydata$conc and mydata$pSurface_cult
## z = 1.8471, p-value = 0.06473
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
## tau
## 0.2531458

La p-value obtenue est très similaire à celle obtenue avec la méthode de Spearman, et se rapproche du seuil de significativité.

Au final compte tenu que les p-values des test de Spearman et de Kendall sont proches de 0.05, que la taille de l’échantillon est faible , et qu’il s’agit des données de terrain, je dirai qu’il semble exister un lien monotone croissant entre les deux variables.

4.2 Présentation des résultats de la corrélation de Spearman

Puisque la corrélation ne repose pas sur un modèle, il n’est pas possible de la représenter graphiquement. Les résultats d’un test de corrélation sont généralement écrits en toutes lettres , ou présentés dans un tableau : on rapporte la valeur du coefficient de corrélation, ainsi que son intervalle de confiance à 95%, et la p-value.

Pour obtenir l’intervalle de confiance à 95% du coefficient de corrélation de Spearman, le plus simple est d’utiliser la fonction slipper_ci du package `slipper`.

 

library(slipper)

set.seed(789) # permet d'obtenir des résultats reproductibles

mydata %>%
slipper_ci(cor.test(conc, pSurface_cult, method="spearman")$estimate)

## ci_low ci_high
## 1 0.005913665 0.6198916

L’intervalle de confiance à 95% du coefficient de corrélation de Spearman est [0.0059 ; 0.62]. 

5 Conclusion

Comme nous l’avons vu ici, lorsque les données s’éloignent du cas idéal, la régression linéaire simple peut poser de nombreuses difficultés. De ce fait, il est souvent plus avantageux, en termes de temps et de robustesse des résultats, d’utiliser l’approche de la corrélation de Spearman. Néanmoins cette méthode n’est utilisable que si la caractérisation de l’évolution d’une des variables en fonction de la seconde, n’est pas au centre des questions.

6 Exprimez vous

Et vous, quelle attitude adoptez face aux pvalues ? Etes vous de ceux qui considèrent qu’une pvalue comprise entre 0.05 et 0.1 témoigne d’une faible preuve, ou de ceux qui refusent de conclure si la pvalue est entre 0.03 et 0.07 ? Utilisez vous une autre règle de décision ? Si vous avez d’autres remarques, n’hésitez pas m’en faire part dans les commentaires, elles sont les bienvenues.

Dans tous les cas, j’espère que ce tutoriel aidera le plus grand nombre à évaluer, avec le logiciel R, si deux variables numériques continues sont liées. Si cet article vous a plu, ou vous a été utile, et si vous le souhaitez, vous pouvez soutenir ce blog en faisant un don sur sa page Tipeee 🙏👉 Cliquez ici pour soutenir le blog Statistiques et Logiciel RCrédit photo : Christian Mayrhofer. Retrouvez la première partie de cet article, ici !