Statistiques Appliquées aux Géosciences

DEUST 1 Géosciences – Travaux Dirigés n°2 (Correction)

Statistiques Bivariées avec R

Objectifs de la séance :

  • Comprendre le concept de statistiques bivariées ;
  • Étudier la relation entre deux variables quantitatives ;
  • Utiliser R pour calculer et interpréter les corrélations, les ajustements statistiques et visualiser les données bivariées ;
  • Analyser les résultats pour interpréter les phénomènes géoscientifiques.

Pré-requis pour la séance :

  • La compréhension totale du TD1 (la partie “statistique” et la partie “programmation R”).

1. Rappels sur l’analyse univariée de données quantitatives

Pour cette partie, vous allez utiliser le jeu de données rock qui contient des mesures de 48 échantillons de roches provenant d’un réservoir pétrolier. Ce jeu de données contient des informations sur la perméabilité des roches, leur surface, leur périmètre et leur forme.

area peri shape perm
4990 2791.90 0.0903296 6.3
7002 3892.60 0.1486220 6.3
7558 3930.66 0.1833120 6.3
  • La colonne area concerne la surface des pores (en nombre de pixels par rapport à la zone composée de \(256 \times 256\) pixels) ;
  • La colonne peri concerne le périmètre (en pixels) ;
  • La colonne shape concerne la forme \((\frac{perimeter}{\sqrt{area}})\) ;
  • La colonne perm concerne la perméabilité (exprimé en millidarcy).

Pour vérifier que vous avez bien assimilé le TD1, les exercices dans cette partie reprennent les mêmes analyses descriptives que durant le TD1, mais avec un jeu de données différent. Bien évidemment, comme pour le TD1, il en convient de charger en amont le jeu de données.

1.1. Combien y a-t-il d’observations dans ce jeu de données ? Pour combien de variables ? À déterminer via une commande R, évidemment.

Indication : Utilisez la fonction dim().

dim(rock)
## [1] 48  4

Le jeu de données rock contient 48 observations décrites par 4 variables.

1.2. Calculer les médiane, moyenne, minimum, maximum et écart-types de surface, périmètre et perméabilité des roches.

c(median(rock$area), mean(rock$area), min(rock$area), max(rock$area), sd(rock$area))
## [1]  7487.000  7187.729  1016.000 12212.000  2683.849
c(median(rock$peri), mean(rock$peri), min(rock$peri), max(rock$peri), sd(rock$peri))
## [1] 2536.195 2682.212  308.642 4864.220 1431.661
c(median(rock$perm), mean(rock$perm), min(rock$perm), max(rock$perm), sd(rock$perm))
## [1]  130.5000  415.4500    6.3000 1300.0000  437.8182

1.3. Affichez un résumé des statistiques descriptives du jeu de données rock.

summary(rock)
##       area            peri            shape              perm        
##  Min.   : 1016   Min.   : 308.6   Min.   :0.09033   Min.   :   6.30  
##  1st Qu.: 5305   1st Qu.:1414.9   1st Qu.:0.16226   1st Qu.:  76.45  
##  Median : 7487   Median :2536.2   Median :0.19886   Median : 130.50  
##  Mean   : 7188   Mean   :2682.2   Mean   :0.21811   Mean   : 415.45  
##  3rd Qu.: 8870   3rd Qu.:3989.5   3rd Qu.:0.26267   3rd Qu.: 777.50  
##  Max.   :12212   Max.   :4864.2   Max.   :0.46413   Max.   :1300.00

1.4. Affichez la distribution de la perméabilité des roches. Commentez.

par(mfrow=c(1,2))

hist(rock$perm, 
     main="Histogramme des perméabilités", # titre de la figure
     ylab="Occurrences", # axes des ordonnées
     xlab=substitute(paste("perméabilité (", italic("mD"), ")")) # axes des abscisses
     )

boxplot(rock$perm, 
        main="Boxplot des perméabilités", 
        ylab=substitute(paste("perméabilités (", italic("mD"), ")"))
        )

À partir de l’histogramme et du boxplot, on observe que sur ce jeu de données, la distribution de perméabilité des roches est fortement déséquilibrée (asymétrique) avec une majorité des roches à moins de 200 mD.

2. Analyse bivariée de données quantitatives

Pour cette partie, vous allez reprendre le jeu de données quakes du TD précédent.

lat long depth mag stations
-20.42 181.62 562 4.8 41
-20.62 181.03 650 4.2 15
-26.00 184.10 42 5.4 43

2.1. À l’aide de la fonction cov(), calculer la covariance entre la magnitude (mag) et la profondeur (depth). Que pouvez-vous en conclure ?

Rappel : La covariance entre deux variables \(X\) et \(Y\) est définie par :
\(\qquad\qquad\)Cov\((X,Y) = \frac{1}{n}\sum\limits_{i=1}^{n} (x_i - \bar{x})(y_i - \bar{y})\).

cov(quakes$mag, quakes$depth) # pareil que l'inverse : cov(quakes$depth, quakes$mag)
## [1] -20.02209

La covariance entre la magnitude et la profondeur est différente de 0, ce qui implique que ces deux variables ne sont pas indépendantes.

Si deux variables \(X\) et \(Y\) sont indépendantes, alors Cov\((X,Y) = 0\).
Attention, la réciproque –si Cov\((X,Y) = 0\), alors \(X\) et \(Y\) sont indépendantes– n’est pas toujours valide.
En revanche, la contraposée est valide : Si Cov\((X,Y) \neq 0\), alors \(X\) et \(Y\) ne sont pas indépendantes.

2.2. Calculer la covariance entre la magnitude d’un séisme et le nombre de stations ayant enregistrées ce séisme. Que pouvez-vous en conclure ?

cov(quakes$mag, quakes$stations) # pareil que l'inverse : cov(quakes$stations, quakes$mag)
## [1] 7.508181

La covariance entre la magnitude d’un séisme et le nombre de stations ayant enregistrées ce séisme est différente de 0, ce qui implique que ces deux variables ne sont pas indépendantes.

2.3. À partir de ces deux covariances calculées, quelles relations pouvez-vous établir concernant la magnitude d’un séisme par rapport à sa profondeur et au nombre de stations ayant enregistrées ce séisme ?

La covariance entre la magnitude et la profondeur est négative, tandis que la covariance entre la magnitude et le nombre de stations ayant enregistrées le séisme est positive. Cela indique que la magnitude a une dépendance négative avec la profondeur tandis qu’elle a une dépendance positive avec le nombre de stations ayant enregistrées le séisme. Nous allons chercher à confirmer cela avec la question suivante…

2.4. À l’aide de la fonction cor(), calculer le coefficient de corrélation de Pearson pour les couples (mag, depth) et (mag, stations). Que pouvez-vous en conclure ?

Rappel : Le coefficient de corrélation de Pearson est défini par :
\(\qquad\qquad\rho(X, Y) = \rho(Y, X) = \frac{\text{Cov}(X,Y)}{\sigma_{X} . \sigma_{Y}}\)
avec \(\bar{x}\) et \(\bar{y}\), la moyenne observée des variables \(x\) et \(y\), respectivement.
Il est borné entre -1 et 1. Il s’agit d’une forme normalisée de la covariance.

La covariance donne seulement une idée de la relation entre deux variables, mais elle dépend fortement des unités de ces deux variables.

Prenons l’exemple du coût des années d’études en Francs et en milliers de Francs avec la durée d’études en années et en mois. Intuitivement, la relation entre coût et durée est la même. Cependant, si l’on considère les différentes unités (francs, années), (francs, mois), (milliers de francs, années) et (milliers de francs, mois), les covariances calculés seront différentes, alors que la relation ne diffère pas.

Il devient donc nécessaire de normaliser l’équation de la covariance afin de pouvoir interpréter l’intensité et la direction de la relation, sans se soucier de l’unité des variables.

cor(quakes$mag, quakes$depth) # pareil que l'inverse : cor(quakes$depth, quakes$mag)
## [1] -0.2306377

La corrélation entre la magnitude et la profondeur est de -0.2306, ce qui indique une faible corrélation négative. Cela signifie que lorsque la magnitude augmente, la profondeur aura tendance à diminuer légèrement, mais elle demeure faible, indiquant qu’il y a beaucoup de variabilité qui n’est pas expliquée par cette relation.

cor(quakes$mag, quakes$stations) # pareil que l'inverse : cor(quakes$stations, quakes$mag)
## [1] 0.8511824

La corrélation entre la magnitude d’un séisme et le nombre de stations ayant enregistrées ce séisme est de 0.8512, ce qui indique une corrélation positive. Un corrélation de 85.12% indique une très forte corrélation entre les deux variables. Ainsi, lorsque la magnitude d’un séisme augmente, le nombre de stations ayant enregistrées ce séisme a tendance à augmenter, et ce de manière quasi proportionnelle.

Attention, une corrélation (même très forte) n’implique pas nécessairement la causalité ! La corrélation mesure l’intensité et la direction d’une relation linéaire entre deux variables pour indiquer si les variables ont tendance à varier ensemble. Quant à elle, la causalité indique qu’une des variables provoque le changement dans l’autre.

Pour faire simple, la corrélation indique ce qui se passe “ensemble”, alors que la causalité indique le changement d’une variable provoqué par une autre (la relation ne va que dans un sens).

2.5. La méthode des moindres carrés ordinaire (MCO)

La méthode des moindres carrés ordinaire (en anglais, Ordinary Least Squares -OLS-) est une technique (parmi d’autres) d’ajustement statistique qui consiste à minimiser la somme des carrés des écarts entre les valeurs observées et les valeurs prédites par un modèle. Elle est souvent utilisée dans les modèles de régression, notamment la régression linéaire.

Principe : Supposons que nous avons un ensemble de données sous la forme de \(n\) paires \((x_i, y_i)\), \(y_i\) est la variable à expliquer en fonction de \(x_i\), la variable explicative. À l’aide de la méthode des moindres carrés ordinaire, nous allons chercher à ajuster une droite de la forme : \(y = ax + b\), où \(a\) est la pente de la droite et \(b\) est l’ordonnée à l’origine (intercept).
Ainsi, la pente et l’ordonnée sont respectivement définies par :

\(\qquad\qquad a = \frac{\text{Cov}(X,Y)}{\sigma^2_X}\qquad\) avec \(\sigma^2_X\) la variance de \(X\)

\(\qquad\qquad b = \bar{Y} - \frac{\text{Cov}(X,Y)}{\sigma^2_X}\bar{X}\quad\) avec \(\bar{X}\) et \(\bar{Y}\), la moyenne observée des variables \(X\) et \(Y\), respectivement.

Pourquoi minimiser la somme des carrés des écarts ? En minimisant la somme des carrés des écarts, la méthode des MCO cherche à rendre les erreurs aussi petites que possible, tout en pénalisant les erreurs les plus importantes (d’où l’élévation au carré). Cela permet d’obtenir une solution unique pour les paramètres \(a\) et \(b\) et de garantir que la droite ajustée est optimale dans le sens des moindres carrés.

Attention, la méthode des moindres carrés ordinaire doit vérifier plusieurs hypothèses au préalable :

  • Linéarité : La relation entre la variable dépendante et les variables indépendantes doit être linéaire.
  • Indépendance des erreurs : Les erreurs doivent être indépendantes les unes des autres (pas d’autocorrélation).
  • Homoscédasticité : La variance des erreurs doit être constante pour toutes les valeurs des variables indépendantes (pas de phénomène d’hétéroscédasticité).
  • Normalité des erreurs : Les erreurs doivent suivre une distribution normale (surtout pour les tests d’hypothèse).

Ici, on supposera que ces hypothèses sont vérifiées.

2.5.1. À l’aide de la fonction lm(), réalisez un modèle de régression linéaire afin d’expliquer la magnitude d’un séisme par le nombre de stations ayant enregistrées ce séisme.

Indication : Écrite telle quelle, la fonction lm(Y ~ X, data=dataset) réalise une régression linéaire afin d’expliquer la variable Y en fonction de la variable explicative X, avec X et Y, deux variables (colonnes) du jeu de données dataset. Pour faire simple, lm(Y ~ X) est la manière sur R de déterminer l’équation \(y = f(x) = ax + b\).

# Créer un modèle de régression linéaire
# On cherche à prédire 'mag' en fonction de (~) 'stations' 
# à partir du jeu de données 'quakes'
modele <- lm(mag ~ stations,
             data = quakes)

# Afficher un résumé du modèle
summary(modele)
## 
## Call:
## lm(formula = mag ~ stations, data = quakes)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.66269 -0.15124 -0.00428  0.13661  0.95227 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.0972676  0.0122067  335.66   <2e-16 ***
## stations    0.0156542  0.0003056   51.23   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2115 on 998 degrees of freedom
## Multiple R-squared:  0.7245, Adjusted R-squared:  0.7242 
## F-statistic:  2625 on 1 and 998 DF,  p-value: < 2.2e-16
plot(x = quakes$stations, # variable en abscisses
     y = quakes$mag, # variable en ordonnée
     main = "Magnitude des séismes en fonction du nombre de stations", # titre
     xlab = "Nombre de stations", ylab = "Magnitude" # label en abscisses
     )
abline(lm(mag ~ stations, data = quakes), # fonction de la droite
       col = "red", # couleur de la droite (visuel)
       lwd=2 # épaisseur de la droite (visuel)
       )

2.5.2. À partir du résumé que vous afficherez, donnez l’équation de la droite. Commentez.

À partir du résumé du modèle, on a pour équation de la droite :
\(\qquad\qquad y = 0.0156542x + 4.0972676\).

L’intercept, égal à 4.0972676, est la valeur estimée de la magnitude d’un séisme lorsque le nombre de stations vaut 0. Bien que le nombre de stations ne puisse pas être égal à 0 dans un contexte réel, cet intercept (qui est la constante) est nécessaire pour l’ajustement de la droite.

Pour chaque station supplémentaire enregistrant un séisme, la magnitude estimée augmente en moyenne de 0.0156542 unité.

2.5.3. Déterminez le coefficient de détermination du modèle, ainsi que sa valeur ajustée. Commentez.

Le coefficient de détermination \((R^2)\) est une mesure indiquant la qualité du modèle de régression. Il est borné entre 0 et 1. Plus \(R^2\) tend vers 1, mieux le modèle est ajusté. En d’autres termes, plus \(R^2\) tend vers 1, mieux l’équation de droite reflète la tendance du jeu de données. Ce coefficient de détermination est défini par :

\(\qquad\qquad R^2 = 1 - \frac{\sum\limits_{i=1}^{n} (y_i - \hat{y_i})^2}{\sum\limits_{i=1}^{n} (y_i - \bar{y})^2}\qquad\) avec

  • \(y_i\) la valeur observée au point \(x_i\) (la valeur que l’on retrouve dans le jeu de données au point \(x_i\)) ;
  • \(\hat{y_i}\) la valeur estimée au point \(x_i\) (la valeur déterminée par l’équation de la droite au point \(x_i\)) ;
  • \(\bar{y}\) la moyenne observée de la variable \(y\).

Il existe un calcul plus “juste” du coefficient de détermination : le coefficient de détermination ajusté \((R^2_{\text{ajusté}})\). Cette variante permet de prendre en compte le nombre de variable et de la taille de l’échantillon (jeu de données). Ce coefficient de détermination ajusté est défini par :

\(\qquad\qquad R^2_{\text{ajusté}} = 1 - \left(\frac{(1-R^2)\times(n-1)}{n-p-1}\right) \qquad\) avec

  • \(n\) la taille de l’échantillon (jeu de données) ;
  • \(p\) le nombre de variables explicatives intégrées dans le modèle.

À partir du résumé du modèle (summary), on observe que le nombre de stations explique une grande partie de la variabilité de la magnitude, puisque son \(R^2\) est égal à 72.45%. En statistiques, le coefficient de détermination est aussi connu sous le nom de la variance expliquée (aussi la variance totale expliquée).