Objectifs

  1. Évaluer et déchiffrer les paramètres d’une analyse de régression linéaire en prenant en compte plusieurs variables catégorielles et numériques, tout en examinant leur impact sur la variable dépendante.

  2. Expliquer le sens profond d’une interaction entre deux variables, tout en élucidant la signification de son coefficient dans le contexte de la modélisation statistique.

  3. Employez le package ‘emmeans’ pour effectuer des comparaisons des moyennes de réponses entre les différents niveaux d’une variable catégorielle, en mettant en lumière les différences significatives et leurs implications.

  4. Démontrer une compréhension approfondie des techniques de normalisation des prédicteurs dans une régression linéaire multiple, en mettant en avant l’importance de cette procédure pour améliorer la stabilité et l’interprétabilité du modèle, tout en minimisant les biais potentiels liés à l’échelle des variables.

I - Régression Linéaire Multiple

Le modèle de régression linéaire multiple représente la relation entre une variable réponse et m prédicteurs \(x-1\), \(x_2\), \(...\) , \(x_m\).

\[y = \beta_0 + \beta_1x_1 + \beta_2x_2 + \ldots + \beta_mx_m + \epsilon = \beta_0 + \sum_{i=1}^m \beta_ix_i + \epsilon\]

Comme dans le cas de la régression linéaire simple, les coefficients \(\beta\) peuvent être calculés à partir de la méthode des moindres carrés. Dans ce modèle, chaque coefficient \(\beta_i\) (sauf \(\beta_0\) ) est la dérivée partielle de \(y\) par rapport à un prédicteur \(x_i\). En d’autres mots, ce coefficient représente la différence moyenne de \(y\) associée à une différence d’une unité de \(x_i\) et aucune différence au niveau des autres prédicteurs.

Un modèle de régression peut inclure plusieurs prédicteurs catégoriels ou numériques. Pour ce cours-ci, nous présenterons des exemples incluant:

  • un prédicteur catégoriel et un prédicteur numérique (dans un contexte expérimental, ce modèle est dénommé analyse de la covariance ou ANCOVA);

  • deux prédicteurs catégoriels (ANOVA à deux facteurs);

  • deux prédicteurs numériques.

Analyse de la covariance

Le tableau de données compensation.csv est tiré du livre de Crawley, Statistics: An introduction using R. Il contient des données sur la masse des graines produites par une espèce de plante (Fruit) en fonction de la taille des racines (Root) et selon que la plante subisse ou non un pâturage (Grazing).

compensation <- read.csv(file.choose())
head(compensation)
##    Root Fruit  Grazing
## 1 6.225 59.77 Ungrazed
## 2 6.487 60.98 Ungrazed
## 3 4.919 14.73 Ungrazed
## 4 5.130 19.28 Ungrazed
## 5 5.417 34.25 Ungrazed
## 6 5.359 35.53 Ungrazed

Visualisation de nuage de points de la base de données concernée

# install.packages("magick")
# install.packages("gganimate")
# install.packages("gifsky")
library(gganimate)
library(ggplot2)
# Créez un nuage de points animé
p <- ggplot(compensation, aes(x = Root, y = Fruit)) +
  geom_point() +
  labs(title = "Nuage de points animé de la base de données Compensation",
       x = "Root", y = "Fruit") +
  transition_states(Grazing, transition_length = 2, state_length = 1)

# Animation du nuage de points
animation <- p + enter_fade() + exit_fade()

anim <- animate(animation, nframes = 50, renderer = gifski_renderer())
knitr::include_graphics(anim)

str(compensation)
## 'data.frame':    40 obs. of  3 variables:
##  $ Root   : num  6.22 6.49 4.92 5.13 5.42 ...
##  $ Fruit  : num  59.8 61 14.7 19.3 34.2 ...
##  $ Grazing: chr  "Ungrazed" "Ungrazed" "Ungrazed" "Ungrazed" ...
summary(compensation)
##       Root            Fruit          Grazing         
##  Min.   : 4.426   Min.   : 14.73   Length:40         
##  1st Qu.: 6.083   1st Qu.: 41.15   Class :character  
##  Median : 7.123   Median : 60.88   Mode  :character  
##  Mean   : 7.181   Mean   : 59.41                     
##  3rd Qu.: 8.510   3rd Qu.: 76.19                     
##  Max.   :10.253   Max.   :116.05

Inspectons d’abord les données.

library(ggplot2)
ggplot(compensation, aes(x = Root, y = Fruit, color = Grazing)) +
  geom_point(size = 3, alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  theme_minimal() +
  labs(
    title = "Dispersion entre la croissance des racines et des fruits",
    x = "Croissance des racines",
    y = "Croissance des fruits",
    color = "Type de Pâturage"
  ) +
  theme(
    plot.title = element_text(size = 20, hjust = 0.5),
    axis.title = element_text(size = 14),
    legend.position = "top"
  )

# Ou :
# ggplot(comp, aes(x = Root, y = Fruit, color = Grazing)) +
#    geom_point() +
#   scale_color_brewer(palette = "Dark2")

Le graphique montre bien l’existence d’une relation linéaire entre la taille des racines et la production de graines, ainsi que l’effet du traitement: pour la même taille des racines, le pâturage réduit la production de graines. Notez que si on n’avait pas mesuré les racines, on pourrait croire que l’effet du pâturage est positif.

ggplot(compensation, aes(x = Grazing, y = Fruit, fill = Grazing)) +
  geom_boxplot(width = 0.6, outlier.shape = NA) +
  scale_fill_brewer(palette = "Dark2") +
  theme_minimal() +
  labs(
    title = "Distribution de la croissance des fruits par type de pâturage",
    x = "Type de Pâturage",
    y = "Croissance des Fruits",
    fill = "Type de Pâturage"
  ) +
  theme(
    plot.title = element_text(size = 20, hjust = 0.5),
    axis.title = element_text(size = 14),
    legend.position = "top"
  )

Cela est dû au fait que les plantes subissant le pâturage avaient (en moyenne) de plus grandes racines au départ. La taille des racines est donc une variable confondante, c’est-à-dire qu’elle est corrélée à la fois avec la variable réponse et avec le traitement dont on cherche à estimer l’effet sur cette réponse. Il faut donc l’inclure dans le modèle pour bien estimer l’effet du pâturage.

Voici un modèle linéaire où l’effet des deux prédicteurs est additif, tel qu’indiqué par le signe + dans la formule du modèle en R:

mod_compensation <- lm(Fruit ~ Grazing + Root, data = compensation)
summary(mod_compensation)
## 
## Call:
## lm(formula = Fruit ~ Grazing + Root, data = compensation)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.1920  -2.8224   0.3223   3.9144  17.3290 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -127.829      9.664  -13.23 1.35e-15 ***
## GrazingUngrazed   36.103      3.357   10.75 6.11e-13 ***
## Root              23.560      1.149   20.51  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.747 on 37 degrees of freedom
## Multiple R-squared:  0.9291, Adjusted R-squared:  0.9252 
## F-statistic: 242.3 on 2 and 37 DF,  p-value: < 2.2e-16

Interprétation des résultats

Si \(x_1\) est la variable de pâturage (0 = Grazed, 1 = Ungrazed selon le codage par défaut dans R) et que \(x_2\) est la taille des racines, l’expression mathématique de ce modèle est:

\(y = \beta_0 + \beta_1x_1 + \beta_2x_2 + \epsilon\)

Pour simplifier l’interprétation des coefficients, on peut séparer le cas avec pâturage (\(x_1\)=0):

\(y = \beta_0 + \beta_2 x_2 + \epsilon\)

du cas sans pâturage (x1=1):

\(y = \beta_0 + \beta_1 + \beta_2 x_2 + \epsilon\)

On peut maintenant intepréter les coefficients comme suit:

  • \(\beta_0\) (Intercept dans le tableau sommaire) est l’ordonnée à l’origine de la droite Fruit vs. Root avec pâturage.

  • \(\beta_1\) (GrazingUngrazed) est l’effet de l’absence de pâturage sur l’ordonnée à l’origine de Fruit vs. Root.

  • \(\beta_2\) (Root) est la pente de la droite Fruit vs. Root avec ou sans pâturage.

Puisque la pente est la même avec ou sans pâturage, le coefficient \(\beta_1\) correspond à une translation sur l’axe des \(y\) de la droite de régression. Ce modèle des effets additifs d’un traitement et d’une variable numérique est donc représenté par deux droites parallèles, ce qui correspond assez bien à notre visualisation des données. En outre, la valeur du \(R^2\) (0.93) indique que le modèle explique une grande partie de la variation observée dans les données.

Même une grande valeur de \(R^2\) ne signifie pas nécessairement que le modèle est approprié. Il faut toujours observer les graphiques de diagnostic. Excepté la présence de quelques valeurs extrêmes dans le diagramme quantile-quantile, les suppositions semblent bien respectées.

Nuage de points du modèle avec prédictions de la variable indépendante

library(ggplot2)
library(gganimate)
library(gifski)

# Créez une copie des données pour éviter de modifier la base de données d'origine
compensation_copy <- compensation

# Créez une nouvelle colonne pour stocker les prédictions du modèle
compensation_copy$Predicted <- predict(mod_compensation)

# Créez un graphique de nuage de points animé
p <- ggplot(compensation_copy, aes(x = Root, y = Fruit, size = Predicted)) +
  geom_point() +
  labs(title = "Nuage de points animé avec prédictions de Fruit",
       x = "Root", y = "Fruit", size = "Predicted") +
  transition_states(Grazing, transition_length = 2, state_length = 1)

# Animation du nuage de points
animation <- p + enter_fade() + exit_fade()

animate(animation, nframes = 50, renderer = gifski_renderer())

Visualisation graphique des hypothèses du modèle

par(mfrow=c(2,2))
plot(mod_compensation)

Notez que le numéro de la rangée du tableau de données est indiqué à côté de certains points extrêmes, pour faciliter l’identification de points problématiques.

Le test F rapporté au bas du sommaire des résultats de lm correspond à l’hypothèse nulle d’absence d’effet pour tous les prédicteurs.

On peut aussi obtenir un tableau d’ANOVA conventionnel en appliquant la fonction anova au résultat de lm.

anova(mod_compensation)
## Analysis of Variance Table
## 
## Response: Fruit
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## Grazing    1  2910.4  2910.4  63.929 1.397e-09 ***
## Root       1 19148.9 19148.9 420.616 < 2.2e-16 ***
## Residuals 37  1684.5    45.5                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Ou sous forme d'un tableau 
# knitr::kable(anova(mod_compensation))

Ordre des prédicteurs

Les fonctions aov et anova dans R traitent les prédicteurs de façon séquentielle, c’est-à-dire que l’effet de chaque prédicteur est calculé par rapport aux résidus du modèle incluant les prédicteurs précédents. Dans notre exemple, la somme des écarts carrés pour le prédicteur Root est basée sur les résidus du modèle incluant seulement Grazing.

C’est ce qu’on appelle une “somme des écarts carrés de Type I” en statistiques. En particulier, cela signifie que le tableau d’ANOVA ne serait pas nécessairement le même en changeant l’ordre des prédicteurs, ex.: Fruit ~ Root + Grazing. D’autres packages en R permettent de réaliser une ANOVA avec des sommes des écarts carrés de Type II et III, mais nous ne verrons pas ces concepts dans ce cours.

Comme nous avons mentionné plus tôt, les coefficients de la régression linéaire multiple estiment l’effet partiel de chaque prédicteur, c’est-à-dire l’effet d’une différence au niveau de ce prédicteur entre deux cas qui ne diffèrent pour aucun autre prédicteur. Pour cette raison, l’ordre des prédicteurs n’influence pas les estimés obtenus avec lm.

Modèle avec interaction

Le modèle précédent suppose que les effets de la taille des racines et du pâturage sur la masse des graines sont additifs: autrement dit, la différence entre les deux traitements de pâturage est la même pour chaque valeur de Root et la pente de Fruit vs. Root est la même pour les cas avec et sans pâturage.

Pour considérer la possibilité que l’effet d’un prédicteur sur la réponse dépende de la valeur d’un autre prédicteur, nous devons spécifier une interaction entre ces deux prédicteurs. Dans la formule d’un modèle en R, l’interaction est indiquée par un symbole de multiplication * entre les variables au lieu du symbole d’addition +.

mod_compensation_inter <- lm(Fruit ~ Grazing * Root, data = compensation)
summary(mod_compensation_inter)
## 
## Call:
## lm(formula = Fruit ~ Grazing * Root, data = compensation)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.3177  -2.8320   0.1247   3.8511  17.1313 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -125.173     12.811  -9.771 1.15e-11 ***
## GrazingUngrazed        30.806     16.842   1.829   0.0757 .  
## Root                   23.240      1.531  15.182  < 2e-16 ***
## GrazingUngrazed:Root    0.756      2.354   0.321   0.7500    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.831 on 36 degrees of freedom
## Multiple R-squared:  0.9293, Adjusted R-squared:  0.9234 
## F-statistic: 157.6 on 3 and 36 DF,  p-value: < 2.2e-16

Si \(x_1\) est la variable de pâturage (0 = Grazed, 1 = Ungrazed selon le codage par défaut dans R) et que \(x_2\) est la taille des racines, l’expression mathématique de ce modèle est:

\(y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_{12} x_1x_2 + \epsilon\)

L’interaction est donc équivalente à l’ajout d’un nouveau prédicteur au modèle, égal au produit des deux variables qui interagissent. Séparons de nouveau en deux équations selon le traitement:

Avec pâturage (\(x_1\)=0):

\(y = \beta_0 + \beta_2 x_2\)

Sans pâturage (\(x_1\)=1):

\(y = (\beta_0 + \beta_1) + (\beta_2 + \beta_{12}) x_2\)

Pour ce modèle avec interaction, l’interprétation des coefficients change un peu:

  • \(\beta_0\) (Intercept dans le tableau sommaire) est l’ordonnée à l’origine de la droite Fruit vs. Root avec pâturage.

  • \(\beta_1\) (GrazingUngrazed) est l’effet de l’absence de pâturage sur l’ordonnée à l’origine de Fruit vs. Root.

  • \(\beta_2\) (Root) est la pente de la droite Fruit vs. Root avec pâturage.

  • \(\beta_12\) (GrazingUngrazed:Root) est l’effet de l’absence de pâturage sur la pente de la droite Fruit vs. Root.

Le modèle avec interaction est donc équivalent à estimer séparément la droite de régression (ordonnée à l’origine et pente) pour chacun des deux traitements.

Comparé au modèle additif, notez que l’effet de l’absence de pâturage (GrazingUngrazed) a maintenant une erreur-type beaucoup plus élevée et une valeur \(p\) plus grande.

summary(mod_compensation)$coefficients
##                   Estimate Std. Error   t value     Pr(>|t|)
## (Intercept)     -127.82936   9.664095 -13.22725 1.349804e-15
## GrazingUngrazed   36.10325   3.357396  10.75335 6.107286e-13
## Root              23.56005   1.148771  20.50892 8.408231e-22
summary(mod_compensation_inter)$coefficients
##                          Estimate Std. Error    t value     Pr(>|t|)
## (Intercept)          -125.1730569  12.811165 -9.7706222 1.150540e-11
## GrazingUngrazed        30.8057049  16.841823  1.8291194 7.567489e-02
## Root                   23.2403732   1.530771 15.1821314 3.173208e-17
## GrazingUngrazed:Root    0.7560338   2.354111  0.3211547 7.499503e-01

Ceci est dû au fait que l’ordonnée à l’origine, correspondant à Root = 0, se situe loin de l’étendue des données (les valeurs de Root sont toutes entre 4 et 11). Donc, un petit changement de pente au milieu du graphique peut mener à un changement important d’ordonnée à l’origine et l’incertitude du coefficient d’interaction (la différence de pente) se répercute aussi sur l’estimation de la différence d’ordonnée à l’origine.

Le tableau d’ANOVA pour un modèle avec interaction inclut la portion de la somme des écarts carrés due à la variation de chaque prédicteur, ainsi que leur interaction. La portion expliquée par l’interaction est basée sur la différence entre les écarts carrés résiduels du modèle sans interaction et ceux du modèle avec interaction. Dans ce cas-ci, l’effet de l’interaction n’est pas significatif, le modèle additif est donc préférable.

anova(mod_compensation_inter)
## Analysis of Variance Table
## 
## Response: Fruit
##              Df  Sum Sq Mean Sq  F value    Pr(>F)    
## Grazing       1  2910.4  2910.4  62.3795 2.262e-09 ***
## Root          1 19148.9 19148.9 410.4201 < 2.2e-16 ***
## Grazing:Root  1     4.8     4.8   0.1031      0.75    
## Residuals    36  1679.6    46.7                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pourquoi l’effet du pâturage (Grazing) est-il significatif dans le tableau d’ANOVA alors que le coefficient GrazingUngrazed du modèle linéaire ne l’est pas? Dans le tableau d’ANOVA, on teste s’il y a une différence significative de la moyenne de Fruit entre les plantes subissant ou non un pâturage, tandis que le coefficient du modèle linéaire réfère à la différence entre les deux traitements lorsque Root est 0 (ordonnée à l’origine). Dans le cas où les variables interagissent, ces deux tests ne sont pas équivalents, car la différence entre les deux droites (avec ou sans pâturage) dépend de la valeur de Root.

ANOVA à deux facteurs

Exemple

Pour illustrer l’ANOVA à deux facteurs, nous utiliserons d’abord le jeu de données growth.csv provenant du manuel Statistics: An Introduction Using R. L’expérience compare le gain de poids (gain) de 48 animaux suivant trois types de régime alimentaire (diet) avec quatre types de suppléments (supplement). Il y a donc 12 groupes (toutes les combinaisons des 3 régimes et 4 suppléments) de 4 individus chacun.

growth <- read.csv(file.choose())
str(growth)
## 'data.frame':    48 obs. of  3 variables:
##  $ supplement: chr  "supergain" "supergain" "supergain" "supergain" ...
##  $ diet      : chr  "wheat" "wheat" "wheat" "wheat" ...
##  $ gain      : num  17.4 16.8 18.1 15.8 17.7 ...
ggplot(growth, aes(x = supplement, y = gain, color = diet)) +
  geom_point(position = position_dodge(width = 0.3), size = 3, alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  theme_minimal() +
  labs(
    title = "Effet des suppléments sur le gain de poids",
    x = "Supplément",
    y = "Gain de poids",
    color = "Régime"
  ) +
  theme(
    plot.title = element_text(size = 20, hjust = 0.5),
    axis.title = element_text(size = 14),
    legend.position = "right"
  )

À première vue, il semble plausible que les effets du régime et du supplément soient additifs, car la différence entre les régimes est semblable d’un supplément à l’autre et la différence entre les suppléments est semblable d’un régime à l’autre. D’ailleurs, le tableau d’ANOVA du modèle avec interaction ne montre pas d’effet significatif de cette interaction:

aov_growth_inter <- aov(gain ~ diet * supplement, data = growth)
summary(aov_growth_inter)
##                 Df Sum Sq Mean Sq F value   Pr(>F)    
## diet             2 287.17  143.59   83.52 3.00e-14 ***
## supplement       3  91.88   30.63   17.82 2.95e-07 ***
## diet:supplement  6   3.41    0.57    0.33    0.917    
## Residuals       36  61.89    1.72                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Notez qu’il est possible d’utiliser la fonction aov ici car nous n’avons que des variables catégorielles et l’échantillon est équilibré (4 réplicats par combinaison de régime et de supplément).

Voici les résultats du modèle additif. Les deux facteurs ont un effet significatif et le régime explique une plus grande portion de la variance du gain de poids (d’après la somme des écarts carrés) que le supplément.

aov_growth_add <- aov(gain ~ diet + supplement, data = growth)
summary(aov_growth_add)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## diet         2 287.17  143.59   92.36 4.20e-16 ***
## supplement   3  91.88   30.63   19.70 3.98e-08 ***
## Residuals   42  65.30    1.55                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Les graphiques de diagnostic n’indiquent pas de problème:

par(mfrow=c(2,2))
plot(aov_growth_add, which=1:2)

D’après le test des étendues de Tukey, les trois régimes ont un effet différent (blé < avoine < orge). Pour les suppléments, agrimore et supersupp ont un effet plus grand que supergain et control.

TukeyHSD(aov_growth_add)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = gain ~ diet + supplement, data = growth)
## 
## $diet
##                   diff       lwr       upr p adj
## oats-barley  -3.092817 -4.163817 -2.021817 0e+00
## wheat-barley -5.990298 -7.061298 -4.919297 0e+00
## wheat-oats   -2.897481 -3.968481 -1.826481 2e-07
## 
## $supplement
##                           diff        lwr        upr     p adj
## control-agrimore    -2.6967005 -4.0583332 -1.3350677 0.0000234
## supergain-agrimore  -3.3814586 -4.7430914 -2.0198259 0.0000003
## supersupp-agrimore  -0.7273521 -2.0889849  0.6342806 0.4888738
## supergain-control   -0.6847581 -2.0463909  0.6768746 0.5400389
## supersupp-control    1.9693484  0.6077156  3.3309811 0.0020484
## supersupp-supergain  2.6541065  1.2924737  4.0157392 0.0000307

Représentation avec les contrastes

Voici les résultats obtenus pour le même modèle avec lm:

lm_growth_add <- lm(gain ~ diet + supplement, data = growth)
summary(lm_growth_add)
## 
## Call:
## lm(formula = gain ~ diet + supplement, data = growth)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.30792 -0.85929 -0.07713  0.92052  2.90615 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          26.1230     0.4408  59.258  < 2e-16 ***
## dietoats             -3.0928     0.4408  -7.016 1.38e-08 ***
## dietwheat            -5.9903     0.4408 -13.589  < 2e-16 ***
## supplementcontrol    -2.6967     0.5090  -5.298 4.03e-06 ***
## supplementsupergain  -3.3815     0.5090  -6.643 4.72e-08 ***
## supplementsupersupp  -0.7274     0.5090  -1.429     0.16    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.247 on 42 degrees of freedom
## Multiple R-squared:  0.8531, Adjusted R-squared:  0.8356 
## F-statistic: 48.76 on 5 and 42 DF,  p-value: < 2.2e-16

Souvenons-nous que par défaut, R utilise un codage de traitement pour représenter les variables catégorielles dans une régression linéaire, où le premier niveau du facteur (en ordre alphabétique) sert de référence. Ici, barley et agrimore sont les niveaux de référence pour le régime et le supplément, respectivement. Nous pouvons donc interpréter chaque coefficient ainsi:

  • l’ordonnée à l’origine est le gain de poids moyen pour les niveaux de référence (orge et agrimore);

  • les coefficients dietoats et dietwheat donnent la différence moyenne de gain entre le régime correspondant (avoine ou blé) et le régime d’orge;

  • les trois derniers coefficients donnent la différence moyenne de gain entre le supplément correspondant et le supplément agrimore.

Le gain de poids moyen pour toute combinaison d’un régime et d’un supplément peut être obtenue en additionnant les coefficients correspondants. Par exemple, le gain moyen pour un régime d’avoine avec le supplément supergain est de: 26.12 (ordonnée à l’origine) - 3.09 (avoine) - 3.38 (supergain) = 19.65.

Tel que vu au dernier cours, nous pouvons modifier les contrastes pour mieux représenter les questions qui nous intéressent. Le code ci-dessous convertit les deux prédicteurs en facteurs, choisit le groupe témoin control comme référence pour supplement et applique un codage d’effet pour diet.

library(dplyr)
growth <- mutate(growth, diet = as.factor(diet),
                 supplement = relevel(as.factor(supplement), ref = "control")) 
contrasts(growth$diet) <- "contr.sum"
colnames(contrasts(growth$diet)) <- c("barley" , "oats")

lm_growth_add <- lm(gain ~ diet + supplement, data = growth)
summary(lm_growth_add)
## 
## Call:
## lm(formula = gain ~ diet + supplement, data = growth)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.30792 -0.85929 -0.07713  0.92052  2.90615 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         20.39861    0.35994  56.673  < 2e-16 ***
## dietbarley           3.02770    0.25451  11.896 4.93e-15 ***
## dietoats            -0.06511    0.25451  -0.256 0.799333    
## supplementagrimore   2.69670    0.50903   5.298 4.03e-06 ***
## supplementsupergain -0.68476    0.50903  -1.345 0.185772    
## supplementsupersupp  1.96935    0.50903   3.869 0.000375 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.247 on 42 degrees of freedom
## Multiple R-squared:  0.8531, Adjusted R-squared:  0.8356 
## F-statistic: 48.76 on 5 and 42 DF,  p-value: < 2.2e-16

Dans ce cas, nous pouvons interpréter les coeffiicents ainsi:

  • l’ordonnée à l’origine est le gain moyen pour le groupe témoin (control), en faisant la moyenne des trois régimes;

  • les coefficients dietbarley et dietoats donnent la différence moyenne de gain des régimes d’orge et d’avoine comparés à la moyenne des trois régimes. La différence moyenne pour le troisième régime (blé) peut être obtenue en prenant l’opposé de la somme des autres effets: -(3.02 - 0.07) = -2.95.

  • les trois derniers coefficients donnent la différence moyenne de gain entre chaque supplément et le groupe témoin.

Modèle avec interaction

Le jeu de données antibiot.csv contient des mesures de prolifération bactérienne (surface couverte en mm2) en fonction de l’humidité (sec, humide) et de la concentration d’antibiotique (faible, modérée, élevée).

# fileEncoding = "UTF-8" permet de lire les accents correctement
antibiot <- read.csv(file.choose(), fileEncoding = "UTF-8")
str(antibiot)
## 'data.frame':    30 obs. of  3 variables:
##  $ Surface      : num  2.1 2.73 1.86 2.36 2.2 ...
##  $ Humidité     : chr  "sec" "sec" "sec" "sec" ...
##  $ Concentration: chr  "faible" "faible" "faible" "faible" ...

Nous devons manuellement spécifier les niveaux du facteur Concentration pour éviter qu’ils ne soient placés en ordre alphabétique.

antibiot$Concentration <- factor(antibiot$Concentration, 
                                 levels = c("faible", "modérée", "élevée"))
levels(antibiot$Concentration)
## [1] "faible"  "modérée" "élevée"

Voici le graphique de ces données. Est-ce qu’un modèle avec des effets additifs de la concentration d’antibiotique et de l’humidité serait approprié ici?

ggplot(antibiot, aes(x = Concentration, y = Surface, color = Humidité)) +
    geom_point(position = position_dodge(width = 0.3)) + 
    scale_color_brewer(palette = "Dark2")

Ici, il y a une interaction claire entre les deux facteurs. Notamment, les conditions humides sont associées à une plus grande surface bactérienne pour les concentrations faible et modérée d’antibiotiques, mais les conditions sèches ont une plus grande surface bactérienne lorsque la concentration est élevée.

Voici le sommaire et les graphiques de diagnostic pour le modèle de la surface bactérienne en fonction de l’interaction entre les deux facteurs.

aov_antibio <- aov(Surface ~ Concentration * Humidité, antibiot)
summary(aov_antibio)
##                        Df Sum Sq Mean Sq F value   Pr(>F)    
## Concentration           2  15.93   7.965    71.5 7.76e-11 ***
## Humidité                1  20.23  20.228   181.6 1.09e-12 ***
## Concentration:Humidité  2  36.40  18.199   163.4 1.05e-14 ***
## Residuals              24   2.67   0.111                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
par(mfrow=c(2,2))
plot(aov_antibio, which=1:2)

L’interaction entre les 3 catégories de concentration et les 2 catégories d’humidité définit 6 groupes, donc il y a 15 comparaisons possibles pour l’interaction, comme le montre le résultat de TukeyHSD.

TukeyHSD(aov_antibio)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Surface ~ Concentration * Humidité, data = antibiot)
## 
## $Concentration
##                      diff        lwr         upr     p adj
## modérée-faible -0.3939894 -0.7667378 -0.02124113 0.0368807
## élevée-faible  -1.7046765 -2.0774249 -1.33192823 0.0000000
## élevée-modérée -1.3106871 -1.6834354 -0.93793878 0.0000000
## 
## $Humidité
##                 diff       lwr       upr p adj
## sec-humide -1.642264 -1.893794 -1.390734     0
## 
## $`Concentration:Humidité`
##                                     diff          lwr        upr     p adj
## modérée:humide-faible:humide -0.82921989 -1.481887432 -0.1765523 0.0073592
## élevée:humide-faible:humide  -4.22827694 -4.880944489 -3.5756094 0.0000000
## faible:sec-faible:humide     -3.61481768 -4.267485222 -2.9621501 0.0000000
## modérée:sec-faible:humide    -3.57357668 -4.226244229 -2.9209091 0.0000000
## élevée:sec-faible:humide     -2.79589383 -3.448561371 -2.1432263 0.0000000
## élevée:humide-modérée:humide -3.39905706 -4.051724600 -2.7463895 0.0000000
## faible:sec-modérée:humide    -2.78559779 -3.438265333 -2.1329302 0.0000000
## modérée:sec-modérée:humide   -2.74435680 -3.397024340 -2.0916893 0.0000000
## élevée:sec-modérée:humide    -1.96667394 -2.619341482 -1.3140064 0.0000000
## faible:sec-élevée:humide      0.61345927 -0.039208277  1.2661268 0.0740073
## modérée:sec-élevée:humide     0.65470026  0.002032716  1.3073678 0.0489732
## élevée:sec-élevée:humide      1.43238312  0.779715574  2.0850507 0.0000070
## modérée:sec-faible:sec        0.04124099 -0.611426550  0.6939085 0.9999549
## élevée:sec-faible:sec         0.81892385  0.166256308  1.4715914 0.0082690
## élevée:sec-modérée:sec        0.77768286  0.125015314  1.4303504 0.0131278

Nous verrons dans la section suivante une méthode facilitant la visualisation de ces comparaisons.

Le modèle linéaire correspondant à cette ANOVA comporte 6 coefficients:

lm_antibio <- lm(Surface ~ Concentration * Humidité, antibiot)
summary(lm_antibio)
## 
## Call:
## lm(formula = Surface ~ Concentration * Humidité, data = antibiot)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.56688 -0.29550  0.04501  0.16490  0.54423 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                        5.8657     0.1493  39.298  < 2e-16 ***
## Concentrationmodérée              -0.8292     0.2111  -3.928 0.000631 ***
## Concentrationélevée               -4.2283     0.2111 -20.031  < 2e-16 ***
## Humiditésec                       -3.6148     0.2111 -17.125 5.87e-15 ***
## Concentrationmodérée:Humiditésec   0.8705     0.2985   2.916 0.007572 ** 
## Concentrationélevée:Humiditésec    5.0472     0.2985  16.907 7.80e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3338 on 24 degrees of freedom
## Multiple R-squared:  0.9645, Adjusted R-squared:  0.9571 
## F-statistic: 130.3 on 5 and 24 DF,  p-value: < 2.2e-16

  • L’ordonnée à l’origine est la surface moyenne pour les niveaux de référence (faible et humide).

  • Les coefficients Concentrationmodérée et Concentrationélevée donnent la différence de surface moyenne due à l’augmentation de concentration de faible à modérée et de faible à élevée, pour le cas humide.

  • Le coefficient Humiditésec donne la différence de surface moyenne entre les cas sec et humide, pour une concentration faible.

  • Finalement, les coefficients liés à l’interaction montrent la différence entre les surfaces moyennes pour les combinaisons “modérée et sec” et “élevée et sec”, comparées aux moyennes prédites par les effets additifs seulement. Autrement dit, la moyenne de la surface bactérienne pour la combinaison “modérée et sec” est égale à: 5.87 (ordonnée à l’origine) - 0.83 (concentration modérée) - 3.61 (sec) + 0.87 (interaction modérée-sec) = 2.30.

Visualisation des effets avec le package emmeans

L’exemple précédent démontre qu’en présence d’une interaction, il est difficile de calculer la moyenne de la réponse pour une combinaison de traitements donnés. Le package emmeans (estimated marginal means) effectue automatiquement le calcul des moyennes pour chaque combinaison de traitements, ainsi que leur intervalle de confiance.

Ci-dessous, nous appliquons la fonction emmeans au résultat du modèle lm_antibio. Le deuxième argument de la fonction spécifie les prédicteurs à considérer: ceux-ci sont indiqués sous forme de formule comme dans la fonction lm, mais sans variable réponse à gauche du ~.

library(emmeans)
em_antibio <- emmeans(lm_antibio, ~ Concentration * Humidité)
em_antibio
##  Concentration Humidité emmean    SE df lower.CL upper.CL
##  faible        humide     5.87 0.149 24     5.56     6.17
##  modérée       humide     5.04 0.149 24     4.73     5.34
##  élevée        humide     1.64 0.149 24     1.33     1.95
##  faible        sec        2.25 0.149 24     1.94     2.56
##  modérée       sec        2.29 0.149 24     1.98     2.60
##  élevée        sec        3.07 0.149 24     2.76     3.38
## 
## Confidence level used: 0.95

La fonction plot appliquée aux résultats d’emmeans illustre les moyennes et leurs intervalles de confiance.

plot(em_antibio)

Il s’agit d’un graphique ggplot2, donc vous pouvez le personnaliser avec les commandes habituelles.

plot(em_antibio) +
    labs(x = "Surface bactérienne moyenne (millimètres carrés)")

Les intervalles de confiance pour chaque moyenne ne nous permettent pas directement de déterminer si deux moyennes diffèrent de façon significative. Pour ce faire, nous spécifions comparisons = TRUE, ce qui ajoute au graphique des flèches de comparaison, basées sur un test de Tukey. Des flèches qui se recoupent sur l’axe de la variable réponse indiquent que les moyennes ne sont pas significativement différentes (à un seuil \(\alpha\)=0.05 , par défaut).

plot(em_antibio, comparisons = TRUE)

Les comparaisons illustrées ici sont les mêmes que celles obtenues précédemment avec le test des étendues de Tukey, mais la visualisation des effets est simplifée. De plus, la fonction TukeyHSD ne s’applique qu’au résultat de la fonction aov, tandis qu’emmeans s’appliquent à tous les modèles de régression que nous verrons dans ce cours.

Lorsqu’un modèle est additif, nous pouvons estimer les moyennes pour un seul facteur. Dans ce cas, l’estimé indiqué correspond à la réponse moyenne pour chaque niveau du facteur, en prenant la moyenne de tous les autres prédicteurs. Dans l’exemple ci-dessous, nous calculons donc la moyenne du gain de poids pour chaque supplément, en faisant la moyenne des trois régimes.

em_growth_supp <- emmeans(lm_growth_add, ~ supplement)
em_growth_supp
##  supplement emmean   SE df lower.CL upper.CL
##  control      20.4 0.36 42     19.7     21.1
##  agrimore     23.1 0.36 42     22.4     23.8
##  supergain    19.7 0.36 42     19.0     20.4
##  supersupp    22.4 0.36 42     21.6     23.1
## 
## Results are averaged over the levels of: diet 
## Confidence level used: 0.95
plot(em_growth_supp)

Régression avec plusieurs prédicteurs numériques

Exemple

Le tableau de données hills du package MASS (inclus par défaut avec R) contient les records de temps (time, en minutes) pour des courses de vélo en Écosse en fonction de la distance horizontale (dist, en milles) et le dénivelé total du parcours (climb, en pieds).

library(MASS)
str(hills)
## 'data.frame':    35 obs. of  3 variables:
##  $ dist : num  2.5 6 6 7.5 8 8 16 6 5 6 ...
##  $ climb: int  650 2500 900 800 3070 2866 7500 800 800 650 ...
##  $ time : num  16.1 48.4 33.6 45.6 62.3 ...

Pour un tableau de données avec plusieurs variables numériques, la fonction plot affiche une matrice de nuages de points pour chaque paire de variables.

plot(hills)

Les temps records semblent dépendre linéairement de la distance et du dénivelé. (La distance et le dénivelé semblent aussi corrélés, nous y reviendrons au prochain cours.) Nous appliquons donc un modèle linéaire à ces données.

mod_hills <- lm(time ~ dist + climb, hills)
par(mfrow=c(2,2))
plot(mod_hills)

Puisque les rangées de ce tableau de données sont identifiées par des noms (rownames dans R), ces noms apparaissent vis-à-vis les valeurs extrêmes dans les graphiques de diagnostic.

D’après ces graphiques, deux parcours (Knock Hill et Bens of Jura) ont un temps record beaucoup plus long qu’attendu (résidu positif important). Ces mêmes points ont aussi une grande influence sur les coefficients de la régression (d’après le quatrième graphique). Dans ce cas, il faudrait vérifier si ces parcours ont des particularités qui expliquent cette forte différence par rapport au modèle.

Graphiques de diagnostic avec lindia

En plus des graphiques de diagnostic obtenus avec plot, il est utile dans le cas d’une régression multiple de visualiser les résidus en fonction de chacun des prédicteurs. La fonction gg_resX du package lindia (pour linear diagnostics) réalise automatiquement ces graphiques à partir de la sortie du modèle.

# install.packages("lindia")
library(lindia)
gg_resX(mod_hills, ncol = 2) # ncol: nombre de colonnes

La présence d’une tendance dans les résidus par rapport à un prédicteur indiquerait un effet non-linéaire possible pour ce prédicteur.

Le package lindia produit aussi d’autres graphiques d diagnostic semblables à ceux obtenus avec plot. Vous pouvez produire tous les graphiques de diagnostic d’un modèle avec la fonction gg_diagnose. Il s’agit de graphiques de type ggplot2, donc vous pouvez les personnaliser avec les fonctions habituelles.

Normalisation des variables

Regardons le sommaire des résultats du modèle:

summary(mod_hills)
## 
## Call:
## lm(formula = time ~ dist + climb, data = hills)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.215  -7.129  -1.186   2.371  65.121 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -8.992039   4.302734  -2.090   0.0447 *  
## dist         6.217956   0.601148  10.343 9.86e-12 ***
## climb        0.011048   0.002051   5.387 6.45e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.68 on 32 degrees of freedom
## Multiple R-squared:  0.9191, Adjusted R-squared:  0.914 
## F-statistic: 181.7 on 2 and 32 DF,  p-value: < 2.2e-16

La valeur des coefficients signifie qu’en moyenne, chaque mille de distance ajoute 6.22 minutes au temps record tandis que chaque pied de dénivelé ajoute 0.01 minute. Puisque les prédicteurs n’ont pas les mêmes unités, la valeur des coefficients n’est pas indicatrice de l’importance de chaque variable. Dans ce cas-ci, dist varie entre 2 et 28 milles tandis que climb varie entre 300 et 7500 pieds.

Aussi, l’ordonnée à l’origine n’a pas vraiment de sens concret, puisqu’un parcours de longueur 0 n’est pas possible.

Afin de comparer l’influence de différents prédicteurs, il peut être utile de les normaliser ceux-ci, c’est-à-dire de transformer chaque valeur en soustrayant la moyenne et en divisant par l’écart-type. Dans R, la fonction scale effectue automatiquement cette transformation.

hills_scl <- hills
hills_scl[, -3] <- scale(hills_scl[, -3]) # on ne normalise pas la réponse
mod_hills_scl <- lm(time ~ dist + climb, data = hills_scl)
summary(mod_hills_scl)
## 
## Call:
## lm(formula = time ~ dist + climb, data = hills_scl)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.215  -7.129  -1.186   2.371  65.121 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   57.876      2.481  23.331  < 2e-16 ***
## dist          34.348      3.321  10.343 9.86e-12 ***
## climb         17.888      3.321   5.387 6.45e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.68 on 32 degrees of freedom
## Multiple R-squared:  0.9191, Adjusted R-squared:  0.914 
## F-statistic: 181.7 on 2 and 32 DF,  p-value: < 2.2e-16

Pour chaque point, la variable normalisée indique l’écart à la moyenne de la variable originale, exprimé en multiple de l’écart-type de la variable originale. Par exemple, dans cette version du modèle, le coefficient dist indique la différence de temps associée à une augmentation d’un écart-type de la distance horizontale. Les coefficients normalisés représentent ainsi l’effet d’une variable relativement aux écarts typiques observés pour cette variable.

Autre avantage de cette représentation, puisque les prédicteurs normalisés prennent une valeur de 0 à leur moyenne, la valeur de l’ordonnée à l’origine de la régression est la moyenne générale de la réponse (ici, le temps record moyen est d’environ 58 minutes).

La normalisation des prédicteurs ne fait que changer l’échelle des effets estimés. La significativité de l’effet de chaque prédicteur et les prédictions du modèle restent les mêmes.

Interaction entre variables continues

Comment interpréter l’interaction entre deux variables continues? Par exemple:

mod_hills_inter <- lm(time ~ dist * climb, hills_scl)
summary(mod_hills_inter)
## 
## Call:
## lm(formula = time ~ dist * climb, data = hills_scl)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -25.994  -4.968  -2.220   2.381  56.115 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   52.304      2.793  18.728  < 2e-16 ***
## dist          32.776      2.965  11.053 2.78e-12 ***
## climb         10.411      3.742   2.782  0.00911 ** 
## dist:climb     8.793      2.745   3.203  0.00314 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.92 on 31 degrees of freedom
## Multiple R-squared:  0.9392, Adjusted R-squared:  0.9333 
## F-statistic: 159.6 on 3 and 31 DF,  p-value: < 2.2e-16

Comme nous avons vu plus tôt, l’équation d’un modèle à deux variables avec interaction est:

\(y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_{12} x_1x_2 + \epsilon\)

On peut ré-écrire cette équation de deux façons:

\(y = \beta_0 + (\beta_1 + \beta_{12}x_2) x_1 + \beta_2 x_2\)

\(y = \beta_0 + \beta_1 x_1 + (\beta_2 + \beta_{12}x_1) x_2\)

  • \(\beta_0\) est la valeur de \(y\) lorsque \(x_1\)=0 et \(x_2\)=0;

  • \(\beta_1\) est l’effet sur \(y\) d’une augmentation d’une unité de \(x_1\) si \(x_2\)=0;

  • \(\beta_2\) est l’effet sur \(y\) d’une augmentation d’une unité de \(x_2\) si \(x_1\)=0;

  • \(\beta_12\) représente à la fois l’augmentation de la pente de la relation \(y\) vs. \(x_1\) lorsque \(x_2\) augmente d’une unité, et l’augmentation de la pente de la relation \(y\) vs. \(x_2\) lorsque \(x_1\) augmente d’une unité.

La normalisation des variables facilite aussi l’interprétation de ces coefficients en présence d’une interaction: par exemple, si chaque prédicteur a une moyenne de zéro, alors \(\beta_1\) est l’effet de \(x_1\) sur y pour un \(x_2\) moyen.

Il peut être utile de visualiser les prédictions du modèle avec interaction. Ci-dessous, nous créons un tableau de données de prédiction avec expand.grid, qui produit toutes les combinaisons de valeurs à partir des vecteurs dist et climb spécifiés. Pour l’illustration des prédictions avec ggplot, nous convertissons climb en variable catégorielle (facteur) pour obtenir des couleurs distinctes sur le graphique.

# Création d'un jeu de données pour les prédictions
hills_nouv <- expand.grid(dist = seq(-2, 2, 0.2), climb = c(-1, 0, 1))
hills_pred <- predict(mod_hills_inter, newdata = hills_nouv, interval = "confidence")
hills_pred <- cbind(hills_nouv, hills_pred)

# Création du graphique avec ggplot2
ggplot(hills_pred, aes(x = dist, y = fit, color = as.factor(climb), fill = as.factor(climb))) +
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.3) +
  geom_line() +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  labs(
    title = "Prédictions avec intervalles de confiance",
    x = "Distance",
    y = "Prédiction",
    color = "Climb",
    fill = "Climb"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 20, hjust = 0.5),
    axis.title = element_text(size = 14),
    legend.position = "right"
  )

Ce graphique illustre bien l’effet d’une interaction positive (coefficient positif de dist:climb): plus l’une des deux variables augmente, plus l’effet de l’autre variable sur la réponse (la pente de la droite) augmente aussi.

Ici, nous avons utilisé le modèle basé sur les prédicteurs normalisés pour réaliser les prédictions; nous aurions pu prendre un modèle basé sur les prédicteurs originaux afin d’obtenir des échelles plus facilement interprétables pour dist et climb.

Résumé

  • Dans une régression linéaire multiple (sans interaction), le coefficient associé à un prédicteur mesure l’effet d’une différence unitaire du prédicteur sur la réponse, si les autres prédicteurs demeurent constants.

  • Une interaction entre deux prédicteurs signifie que l’effet d’un prédicteur sur la réponse (i.e. la pente de la droite de régression) dépend de la valeur d’un autre prédicteur.

  • Le package emmeans permet d’effectuer des comparaisons multiples pour l’effet d’une variable catégorielle sur une réponse, comme le test des étendues de Tukey, mais applicable à tout modèle de régression.

  • La normalisation des prédicteurs d’une régression (soustraire la moyenne et diviser par l’écart-type) facilite la comparaison des coefficients et l’interprétation de l’ordonnée à l’origine (qui représente la moyenne générale de la variable réponse).

II - Visualisation animée du taux d’inflation durant la période 1982 - 2006

# install.packages("gganimate")
# install.packages("gifski")
library(ggplot2)
library(gganimate)
library(dplyr)
library(gifski)

# Importation des données :
inflationhk <- read.csv2(file.choose()); head(inflationhk)
##   Year Composite.CPI Inflation
## 1 1982          23.8      10.9
## 2 1983          26.1      10.0
## 3 1984          28.4       8.6
## 4 1985          29.4       3.5
## 5 1986          30.4       3.8
## 6 1987          32.1       5.7
# Filtrer les données pour la plage de Year de 1982 à 2006
inflationhk_filtered <- inflationhk %>% filter(Year >= 1982 & Year <= 2006)

# Création du graphique avec ggplot2 en utilisant les données filtrées
anim <- inflationhk_filtered %>% ggplot(aes(Year, Inflation)) + 
  geom_line(color = "darkorange") + geom_point(color = "darkorange") + 
  ggtitle("HK Inflation Rate 1982 - 2006") + ylab("Inflation Rate %") +
  theme(legend.position = "none", panel.background = element_rect(fill = "lightblue")) +
  scale_x_continuous(n.breaks = 10) + scale_y_continuous(n.breaks = 10) +
  transition_reveal(Year) + ease_aes('linear')

animate(anim, renderer = gifski_renderer()) 

anim_save("inflation.gif")

III - Visualisation OLS (Assumptions de Régression Linéaire)

En général, il est suffisant de valider vos hypothèses de modèle de manière programmée, mais parfois il est utile de les visualiser. Voici quelques moyens rapides pour le faire.

Linéarité

La linéarité est probablement l’hypothèse la plus facile à visualiser, car vous pouvez simplement utiliser le code suivant pour créer rapidement un graphique de dispersion.

# Génération de données aléatoires pour x et y
set.seed(123) # Pour la reproductibilité
x <- rnorm(100) # Génère 100 valeurs aléatoires pour x
y <- rnorm(100) # Génère 100 valeurs aléatoires pour y

# Tracer le graphique
plot(x, y, main = "Graphique de valeurs aléatoires", xlab = "X", ylab = "Y")

De plus, on peut modifier l’apparence de vos points en utilisant les options “pch”, “cex” et “col”. “PCH” signifie “Plot Character” et permet d’ajuster le symbole utilisé pour vos points. Les formes de points disponibles sont répertoriées dans l’image ci-dessous.

ggpubr::show_point_shapes()
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

L’option “cex” nous permet d’ajuster la taille du symbole. La valeur par défaut est de 1. Si nous modifions la valeur à 0,75, par exemple, le symbole du graphique sera réduit à 3/4 de la taille par défaut. L’option “col” nous permet d’ajuster la couleur de nos symboles de graphique.

plot(x, y, col=rgb(0.4,0.4,0.8,0.6), pch=16, cex=1.2)

Nous pouvons ajuster les axes avec les options “xlab”, “ylab”, “xaxt” et “yaxt” (parmi d’autres options disponibles). Dans l’exemple suivant, nous allons supprimer complètement les axes.

plot(x, y, col=rgb(0.4,0.4,0.8,0.6), pch=16, cex=1.2, xlab="", ylab="", xaxt="n", yaxt="n")

Enfin, nous pouvons ajouter une ligne de tendance en créant un modèle et en ajoutant les valeurs ajustées au graphique. Nous ajusterons également l’épaisseur de la ligne et la couleur avec les paramètres “lwd” et “col”, respectivement. En utilisant ces deux commandes :

model <- lm(y ~ x) lines(model$fitted.values, col=2, lwd=2)

En alternative, vous pouvez enrichir vos données avec des limites en utilisant la fonction “predict” comme indiqué ci-dessous.

# create your model
model <- lm(y ~ x)

# predict your model
predict_model <- predict(model, interval="predict")
## Warning in predict.lm(model, interval = "predict"): les prédictions sur les données actuelles se réfèrent à des réponses _futures_
# plot your raw data
plot(x, y, col=rgb(0.4,0.4,0.8,0.6), pch=16, cex=1.2, xlab="", ylab="", xaxt="n", yaxt="n")

# get the index of your data
ix <- sort(x, index.return=T)$ix

# add your trendline
lines(x[ix], predict_model[ix, 1], col=2, lwd=2)

# add a shape to represent your upper and lower limits
polygon(c(rev(x[ix]), x[ix]), c(rev(predict_model[ix, 3]), predict_model[ix, 2]), col = rgb(0.7,0.7,0.7,0.4), border = NA)

Multicolinéarité

La première manière de visualiser la multicolinéarité est à travers une matrice de graphiques grâce à la fonction “pairs”. Nous pouvons tester cela avec l’ensemble de données “mtcars” de la manière suivante :

pairs(mtcars, pch=20, lower.panel=NULL, xaxt="n", yaxt="n", col="#FC4E07")

La deuxième manière de visualiser la multicolinéarité est à travers un graphique de corrélation. Tout d’abord, installez la bibliothèque “corrplot,” puis utilisez la fonction “corrplot.”

library(corrplot)
## corrplot 0.92 loaded
corrplot(cor(mtcars), method = "number")

Autocorrélation

Pour visualiser l’autocorrélation, nous pouvons créer un graphique d’autocorrélation à l’aide de la fonction “acf” de la bibliothèque “stats”.

library(stats)

model <- lm(mpg~drat, data=mtcars)
acf(model$residuals, type="correlation")

IV - Exploration de Distribution d’une base de données avec Ridgeline Plots

Un graphique de type “ridgeline,” également connu sous le nom de “joyplot,” est une technique de visualisation qui nous permet d’étudier la distribution d’une variable numérique à travers plusieurs groupes. Dans cet article, nous allons explorer comment créer des graphiques de type “ridgeline” en utilisant R et la bibliothèque ggridges. Nous fournirons des exemples détaillés pour vous aider à comprendre et à utiliser efficacement cet outil de visualisation puissant.

Exemple : Prix des diamants par qualité

Commençons par un exemple simple pour illustrer les bases de la création d’un graphique de type “ridgeline.” Dans ce scénario, nous souhaitons examiner la distribution des prix des diamants en fonction de leur qualité. Voici comment nous pouvons générer le graphique de type “ridgeline” en utilisant la bibliothèque ggridges :

# Load required libraries
# install.packages("ggridges")
library(ggridges)
library(ggplot2)

# Create a ridgeline plot
ggplot(diamonds, aes(x = price, y = cut, fill = cut)) +
  geom_density_ridges(scale = 3, rel_min_height = 0.01, alpha = 0.8, color = "white") +
  theme_ridges() + 
  theme(legend.position = "none")
## Picking joint bandwidth of 458

Dans le code ci-dessus, nous utilisons la fonction ggplot() pour initialiser le graphique et geom_density_ridges() de la bibliothèque ggridges pour créer le graphique de type “ridgeline.” En définissant les esthétiques x et y sur les colonnes price et cut, respectivement, nous pouvons visualiser la relation entre ces variables. L’esthétique fill est positionnée à cut pour représenter différentes catégories avec des couleurs distinctes.

ggplot2movies

Le package ggplot2movies fournit des données de l’IMDB sur un grand nombre de films, y compris leurs durées, dans un tibble nommé “movies”.

# install.packages("ggplot2movies")
library(ggplot2movies)
dim(movies)
## [1] 58788    24
head(movies)
## # A tibble: 6 × 24
##   title      year length budget rating votes    r1    r2    r3    r4    r5    r6
##   <chr>     <int>  <int>  <int>  <dbl> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 $          1971    121     NA    6.4   348   4.5   4.5   4.5   4.5  14.5  24.5
## 2 $1000 a …  1939     71     NA    6      20   0    14.5   4.5  24.5  14.5  14.5
## 3 $21 a Da…  1941      7     NA    8.2     5   0     0     0     0     0    24.5
## 4 $40,000    1996     70     NA    8.2     6  14.5   0     0     0     0     0  
## 5 $50,000 …  1975     71     NA    3.4    17  24.5   4.5   0    14.5  14.5   4.5
## 6 $pent      2000     91     NA    4.3    45   4.5   4.5   4.5  14.5  14.5  14.5
## # ℹ 12 more variables: r7 <dbl>, r8 <dbl>, r9 <dbl>, r10 <dbl>, mpaa <chr>,
## #   Action <int>, Animation <int>, Comedy <int>, Drama <int>,
## #   Documentary <int>, Romance <int>, Short <int>

Un graphique de type “ridgeline” des durées des films pour chaque année :

library(dplyr)
mv12 <- filter(movies, year > 1912)
ggplot(mv12, aes(x = length, y = year, group = year)) +
    geom_density_ridges(scale = 10, size = 0.25, rel_min_height = 0.03) +
    scale_x_continuous(limits = c(0, 200)) +
    scale_y_reverse(breaks = c(2000, 1980, 1960, 1940, 1920)) +
    theme_minimal()
## Picking joint bandwidth of 6.68
## Warning: Removed 250 rows containing non-finite values
## (`stat_density_ridges()`).

Cela montre que depuis le début des années 1960, la durée des longs métrages s’est stabilisée autour d’une distribution centrée autour de 90 minutes.

V - Density Plots

Bases de Density Plots

Les graphiques de densité peuvent être considérés comme des représentations graphiques des histogrammes lissés.

La lissage est contrôlé par un paramètre de largeur de bande qui est analogue à la largeur des bacs de l’histogramme.

La plupart des graphiques de densité utilisent une estimation de la densité par noyau, mais il existe d’autres stratégies possibles ; qualitativement, la stratégie particulière a rarement de l’importance.

Voici un exemple de graphique de densité de la variable “durée de l’eruption de geyser” avec une largeur de bande par défaut :

ggplot(geyser) +
    geom_density(aes(x = duration))

L’utilisation d’une largeur de bande plus petite montre l’accumulation à 2 et 4 minutes :

ggplot(geyser) +
    geom_density(aes(x = duration), bw = 0.05)

Pour un nombre modéré d’observations, une addition utile est un “jittered rug plot” :

ggplot(geyser) +
    geom_density(aes(x = duration)) +
    geom_rug(aes(x = duration, y = 0),
             position =
                 position_jitter(height = 0))

Scalability (Scalabilité Visuelle)

La scalabilité visuelle est très bonne.

Pour la variable de prix des diamants dans les données :

ggplot(diamonds) +
    geom_density(aes(x = price))

Les estimations de densité sont généralement calculées sur une grille de points et interpolées.

Les valeurs par défaut dans R varient de 50 à 512 points.

L’effort de calcul pour une estimation de densité en un point est proportionnel au nombre d’observations.

L’espace de stockage nécessaire pour une image est proportionnel au nombre de points où la densité est estimée.

Multiple Groups (Multiple Groupes)

Pour les données sur les diamants, un graphique de densité des prix facettes en fonction de la qualité montre les distributions conditionnelles des prix aux différents niveaux de qualité :

ggplot(diamonds) +
    geom_density(aes(x = price)) +
    facet_wrap(~ cut)

En associant l’esthétique y à after_stat(count), on affiche la distribution conjointe du prix et de la qualité :

ggplot(diamonds) +
    geom_density(aes(x = price,
                     y = after_stat(count))) +
    facet_wrap(~ cut)

Un graphique de densité empilé est parfois utile, mais il est souvent difficile à lire :

ggplot(diamonds) +
    geom_density(aes(x = price,
                     y = after_stat(count),
                     fill = cut),
                 position = "stack")

Une option intermédiaire : un graphique avec facettes sur l’échelle des comptages, avec un graphique atténué pour l’ensemble des données afin de permettre d’évaluer les proportions de l’ensemble :

ggplot(diamonds) +
    geom_density(aes(x = price, y = after_stat(count)),
                 fill = "lightgrey", color = NA,
                 data = mutate(diamonds, cut = NULL)) +
    geom_density(aes(x = price,
                     y = after_stat(count),
                     fill = cut),
                 position = "stack", color = NA) +
    facet_wrap(~ cut) +
    scale_fill_viridis_d(guide = "none")

Un graphique de densité rempli offre une vue de la distribution conditionnelle de la qualité à différents niveaux de prix :

ggplot(diamonds) +
    geom_density(aes(x = price, y = after_stat(count), fill = cut),
                 position = "fill") +
    ylab(NULL)

Cela s’appelle un graphique de densité conditionnelle, ou un “CD plot”.

Boxplots

Les boîtes à moustaches, ou diagrammes en boîte, fournissent une représentation schématique d’une distribution.

Ils conviennent très bien pour montrer les distributions de plusieurs groupes.

Il existe de nombreuses variations des boîtes à moustaches :

  • La plupart commencent par une boîte du premier au troisième quartile et sont divisées par la médiane.

  • La forme la plus simple ajoute ensuite une moustache du premier quartile au minimum et du troisième quartile au maximum.

  • Plus courant est de dessiner la moustache supérieure jusqu’au point le plus élevé en dessous du troisième quartile +1,5∗IQR, et la moustache inférieure de manière analogue.

  • Les valeurs aberrantes qui se trouvent en dehors de la plage des moustaches sont ensuite représentées directement :

# install.packages("gapminder")
library(gapminder)
library(ggplot2)
ggplot(gapminder) +
    geom_boxplot(aes(x = continent, y = gdpPercap)) +
    xlab(NULL)

Il existe des variantes qui distinguent entre les valeurs aberrantes légères et les valeurs aberrantes extrêmes.

Une variante courante consiste à montrer une estimation approximative de l’intervalle de confiance à 95 % pour la médiane de la population sous forme d’une encoche :

ggplot(gapminder) +
    geom_boxplot(aes(x = continent, y = gdpPercap),
                 notch = TRUE) +
    xlab(NULL)

Une autre variante consiste à utiliser une largeur proportionnelle à la racine carrée de la taille de l’échantillon pour refléter la force des preuves dans les données :

ggplot(gapminder) +
    geom_boxplot(aes(x = continent, y = gdpPercap),
                 notch = TRUE, varwidth = TRUE) +
    xlab(NULL)

Avec des tailles d’échantillon modérées, il peut être utile de superposer les données d’origine, peut-être avec des décalages et un mélange alpha.

Les valeurs aberrantes dans le diagramme en boîte peuvent être désactivées avec l’argument outlier.color = NA pour qu’elles ne soient pas affichées deux fois :

p <- ggplot(gapminder) +
    geom_boxplot(aes(x = continent, y = gdpPercap),
                 notch = TRUE, varwidth = TRUE,
                 outlier.color = NA) +
    xlab(NULL)
p + geom_point(aes(x = continent, y = gdpPercap),
               position =
                   position_jitter(width = 0.1),
               alpha = 0.1)

VI - Violin Plots

Une variante du diagramme en boîte est le graphique en violon :

Le graphique en violon utilise des estimations de densité pour montrer les distributions :

ggplot(gapminder) +
    geom_violin(aes(x = continent, y = gdpPercap)) +
    xlab(NULL)

Par défaut, les “violons” sont mis à l’échelle de manière à avoir la même aire.

Ils peuvent également être mis à l’échelle pour avoir la même hauteur maximale ou pour avoir des aires proportionnelles aux tailles des échantillons.

Cela se fait en ajoutant

  • scale = “width” ou
  • scale = “count”

à l’appel de geom_violin.

Une comparaison entre les boîtes à moustaches et les graphiques en violon :

ggplot(gapminder) +
    geom_boxplot(aes(x = continent, y = gdpPercap)) +
    geom_violin(aes(x = continent, y = gdpPercap),
                fill = NA, scale = "width",
                linetype = 2) +
    xlab(NULL)

Une combinaison de boîtes à moustaches et de graphiques en violon :

ggplot(gapminder) +
    geom_violin(aes(x = continent, y = gdpPercap),
                scale = "width") +
    geom_boxplot(aes(x = continent, y = gdpPercap),
                 width = .1) +
    xlab(NULL)

Il existe d’autres variations, par exemple les graphiques en forme de vase.

Les boîtes à moustaches ne reflètent pas la forme d’une distribution.

Pour les éruptions dans l’ensemble de données “faithful” :

ggplot(faithful) +
    geom_boxplot(aes(y = eruptions, x = "Box")) +
    geom_violin(aes(y = eruptions, x = "Violin")) +
    xlab(NULL)

VII - Swarm Plots

Les graphiques en essaim (swarm plots) montrent l’ensemble des données sous une forme qui permet également de voir la densité.

Il existe plusieurs variations et noms, notamment les graphiques en essaim d’abeilles, les graphiques en essaim de violons, les graphiques en essaim de violons, et les graphiques Sina.

Les graphiques Sina sont disponibles sous la forme de geom_sina dans le package ggforce :

# install.packages("ggforce")
library(ggforce)
ggplot(gapminder,
       aes(x = continent, y = gdpPercap)) +
    geom_sina(size = 0.2) +
    xlab(NULL)

Combiné avec un graphique en violon mis à l’échelle en largeur :

ggplot(gapminder,
       aes(x = continent, y = gdpPercap)) +
    geom_violin(scale = "width") +
    geom_sina(color = "blue",
              size = 0.4,
              scale = FALSE) +
    xlab(NULL)

Efficacité et Scalabilité

  • Les boîtes à moustaches (boxplots) sont très simples et faciles à comparer.

  • Les boîtes à moustaches mettent fortement l’accent sur la moitié centrale des données.

  • Les boîtes à moustaches peuvent ne pas être faciles à comprendre pour un spectateur non averti.

  • Les boîtes à moustaches se redimensionnent assez bien visuellement et computationnellement en fonction du nombre d’observations ; le superposition/stockage des valeurs aberrantes devient un problème pour les ensembles de données plus volumineux.

  • Les graphiques en violon (violin plots) se redimensionnent bien à la fois visuellement et computationnellement en fonction du nombre d’observations.

# install.packages("patchwork")
library(patchwork)
p1 <- ggplot(diamonds) +
    geom_boxplot(aes(x = cut, y = price)) +
    xlab(NULL)
p2 <- ggplot(diamonds) +
    geom_violin(aes(x = cut, y = price)) +
    xlab(NULL)
p1 + p2

  • La scalabilité en termes du nombre de cas pour les graphiques en essaim (swarm plots) ou les graphiques Sina est plus limitée.

  • Le nombre de groupes qui peuvent être gérés pour la comparaison à l’aide de ces graphiques se situe dans la plage de quelques dizaines.

library(lattice)
p1 <- ggplot(barley) +
    geom_boxplot(aes(x = site, y = yield, fill = year)) +
    xlab(NULL)
p2 <- ggplot(barley) +
    geom_violin(aes(x = site, y = yield, fill = year)) +
    xlab(NULL)
p1 + p2

Les axes peuvent être inversés pour éviter la superposition des étiquettes :

library(lattice)
p3 <- p1 + coord_flip() + guides(fill = "none")
p4 <- p2 + coord_flip()
p3 + p4

La création de facettes peut également être utilisée pour organiser des groupes de boîtes à moustaches ou de graphiques en violon.

Par exemple, pour l’espérance de vie par continent au fil des années dans les données de gapminder :

library(dplyr)
ggplot(filter(gapminder,
              year %% 10 == 2,
              continent != "Oceania")) +
    geom_boxplot(aes(x = lifeExp, y = factor(year))) +
    facet_wrap(~ continent, ncol = 1) +
    theme_minimal() +
    theme(text = element_text(size = 12)) +
    theme(strip.text.x = element_text(hjust = 0)) +
    ylab(NULL)

Veuillez trouver ci-après mon adresse électronique Gmail si tant est que vous souhaiteriez me contacter pour une quelconque demande d’information ou si vous avez des questions qui vous taraudent :

Merci!