Plan de match

1. Modèles statistiques
2. Un truc pour ajouter les points derrière des graphs
3. Annotations
4. Faire pivoter les données
5. Explications du travail final

Les données du jour

library(tidyverse)
air <- airquality

Modèle de régression linéaire: Rappels

Démo au tableau.

Modèle de régression linéaire

Supposons que nous sommes intéressés par le modèle suivant:

modele <- lm(Temp ~ Wind + Solar.R + as.factor(Month), data = air)

Explorer les résultats

summary(modele)

## 
## Call:
## lm(formula = Temp ~ Wind + Solar.R + as.factor(Month), data = air)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.7987  -4.0036  -0.6883   4.0873  15.7654 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       70.123209   2.295297  30.551  < 2e-16 ***
## Wind              -0.723387   0.144566  -5.004 1.67e-06 ***
## Solar.R            0.023372   0.005513   4.239 4.06e-05 ***
## as.factor(Month)6 11.958897   1.568042   7.627 3.44e-12 ***
## as.factor(Month)7 15.188732   1.602650   9.477  < 2e-16 ***
## as.factor(Month)8 16.164466   1.626746   9.937  < 2e-16 ***
## as.factor(Month)9 10.227537   1.571535   6.508 1.28e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.856 on 139 degrees of freedom
##   (7 observations deleted due to missingness)
## Multiple R-squared:  0.6136, Adjusted R-squared:  0.597 
## F-statistic: 36.79 on 6 and 139 DF,  p-value: < 2.2e-16

La fonction summary() nous montre les résultats de régression sous forme de tableau, avec des colonnes pour les coefficients, les erreurs type, les valeurs t et les statistiques p. 

L’objet modele

L’objet modele n’est pas organisé sous forme de tableau.

View(modele)

Organiser l’information dans un format tidy

• chaque ligne est une observation
• chaque colonne est une variable
• chaque cellule contient une seule valeur

La fonction tidy()

• transforme les résultats de régression sous forme de tableau
• ligne pour les paramètres (variables), colonnes pour les statistiques
• ajout des intervalles de confiance (conf.int = TRUE)
• intervalles à 95%.

modele_tab <- tidy(modele, conf.int = TRUE)
modele_tab

## # A tibble: 7 × 7
##   term              estimate std.error statistic  p.value conf.low conf.high
##   <chr>                <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
## 1 (Intercept)        70.1      2.30        30.6  1.55e-63  65.6      74.7   
## 2 Wind               -0.723    0.145       -5.00 1.67e- 6  -1.01     -0.438 
## 3 Solar.R             0.0234   0.00551      4.24 4.06e- 5   0.0125    0.0343
## 4 as.factor(Month)6  12.0      1.57         7.63 3.44e-12   8.86     15.1   
## 5 as.factor(Month)7  15.2      1.60         9.48 9.62e-17  12.0      18.4   
## 6 as.factor(Month)8  16.2      1.63         9.94 6.56e-18  12.9      19.4   
## 7 as.factor(Month)9  10.2      1.57         6.51 1.28e- 9   7.12     13.3

Visualiser

• Une fois les informations de notre modèle “nettoyées”, nous pouvons nous en servir pour présenter les résultats de régression sous forme de graphique.

graph_mod <- ggplot(aes(y = term,
                        x = estimate), 
                    data = modele_tab) +
  geom_point() +
  theme_minimal()

graph_mod

Enlever des coefficients

• La constante (intercept) ne m’intéresse pas.
• Vu que nos résultats de régression sont organisés sous forme de tableau, il est plutôt simple de faire des modifications.

modele_tab2 <- modele_tab  %>% 
  filter(!str_detect(term, "Intercept")) %>% 
  mutate(term = case_when(str_detect(term, "Solar")~"Solaire (rad.)",
                          str_detect(term, "Wind")~"Vent (vitesse)",
                          str_detect(term, "9")~"Sept.",
                          str_detect(term, "8")~"Août",
                          str_detect(term, "7")~"Juil.",
                          str_detect(term, "6")~"Juin",
                          TRUE~term))

Modifier notre graphique

• J’aimerais que les noms de variables soient à l’horizontal (sur l’axe des y);
• J’aimerais avoir un point pour le coefficient estimé par le modèle, et une ligne qui traverse le point pour montrer l’intervalle de confiance;
• J’aimerais avoir une ligne verticale à 0, pour montrer quels estimés ne sont pas statistiquement significatifs.

Modifier notre graphique

graph_mod <- ggplot(aes(y = term, x = estimate,
                        xmin = conf.low, xmax = conf.high),
                    data = modele_tab2) +
  geom_pointrange() +
  geom_vline(xintercept = 0) +
  labs(title = "Effet du vent sur la température",
       y = "",
       x = "Coefficient et intervalles de confiance (95%)") +
  theme_minimal()

graph_mod

Visualiser des prédictions

• Il est possible de créer des prédictions pour des observations fictives à partir de notre modèle de régression.
• Voir: http://www.sthda.com/english/articles/40-regression-analysis/166-predict-in-r-model-predictions-and-confidence-intervals/
• Voir le chapitre 6 de Healy: https://socviz.co/modeling.html#modeling

Visualiser des prédictions: créer des observations fictives

• Nous allons créer une fausse banque de données avec 2500 observations (vs 153 pour la banque de données air).
• Nous générons 500 valeurs situées entre le minimum (1.7) et le maximum (20.7) de la variable Wind.
• Ces 500 valeurs sont associées à chacun des 5 mois du modèle (5 à 9, ou mai, juin, juillet, août, septembre), ce qui crée les 2500 observations.
• Toutes les observations ont la même valeur pour la variable Solar.R (fixée à la moyenne).

Visualiser des prédictions: créer des observations fictives

• Il serait possible, bien sûr, de créer des fausses observations plus grandes ou plus petites que la vraie étendue de la variable Wind (1.7 à 20.7).
• Or, nos prédictions ne seraient alors pas réalistes.
• Il est donc préférable de s’en tenir aux valeurs “possibles”.

Visualiser des prédictions: créer des observations fictives

air_pred <- expand.grid(Wind = (seq(from = min(air$Wind),
                                   to = max(air$Wind),
                                   length.out = 500)),
                        Solar.R = mean(air$Solar.R, na.rm = T),
                        Month = c(5:9))

head(air_pred)

##       Wind  Solar.R Month
## 1 1.700000 185.9315     5
## 2 1.738076 185.9315     5
## 3 1.776152 185.9315     5
## 4 1.814228 185.9315     5
## 5 1.852305 185.9315     5
## 6 1.890381 185.9315     5

Visualiser des prédictions: Prédire des résultats pour nos observations fictives

• On utilise ensuite la fonction predict() pour trouver quelle valeur de Y (Temp) est prédite (fit) par notre modèle pour chaque valeur de la variabe X (Wind).
• J’en profite pour convertir le résultat en tableau de données avec as.data.frame.

predictions <- predict(object = modele, 
                       newdata = air_pred, 
                       interval = "predict") %>% 
  as.data.frame()

Visualiser des prédictions: Prédire des résultats pour nos observations fictives

• predict() calcule aussi l’intervalle de confiance autour de chaque prédiction (lwr et upr)
• Nous allons joindre ces prédictions au tableau de données fictives créé plus tôt.

Visualiser des prédictions: Prédire des résultats pour nos observations fictives

prediction_viz <- cbind(air_pred, predictions)

head(prediction_viz)

##       Wind  Solar.R Month      fit      lwr      upr
## 1 1.700000 185.9315     5 73.23913 61.10554 85.37272
## 2 1.738076 185.9315     5 73.21158 61.08056 85.34261
## 3 1.776152 185.9315     5 73.18404 61.05556 85.31252
## 4 1.814228 185.9315     5 73.15650 61.03056 85.28243
## 5 1.852305 185.9315     5 73.12895 61.00555 85.25236
## 6 1.890381 185.9315     5 73.10141 60.98053 85.22229

Visualiser des prédictions: Faire le graph

• Nous visualisons les températures prédites par notre modèle pour différentes valeurs de la vitesse du vent.
• Les prédictions ont été calculées en gardant le niveau de radiation solaire constant (à sa moyenne).
• Les prédictions sont calculées pour chaque mois.

Visualiser des prédictions: Faire le graph

graph_pred <- ggplot(aes(x = Wind, y = fit, 
                         ymin = lwr, ymax = upr), 
                     data = prediction_viz) +
# lignes de prédiction (5 couleurs)
  geom_line(aes(color = as.factor(Month))) +
# geom_ribbon ajoute l'intervalle de confiance 
# (une par mois, c'est pourquoi on précise "group")
  geom_ribbon(aes(group = as.factor(Month)), alpha = .1) +
  scale_color_discrete(name = "Mois", 
                      labels = c("Mai", "Juin", "Juillet", 
                                 "Août", "Septembre")) +
  labs(title = "Température prédite en fonction du vent",
       y = "Température prédite\n(avec intervalle de 95%)",
       x = "Vitesse du vent") +
  theme_minimal()

Visualiser des prédictions: Faire le graph

graph_pred

Visualiser des prédictions: Montrer les vraies observations

• Nous pourrions aussi faire apparaître les vraies observations derrière les lignes de prédiction.
• Pour ce faire, j’ajoute un étage geom_point qui utilise les données de la banque de données originales, et non la banque de données fictives.

Visualiser des prédictions: Montrer les vraies observations

graph_pred2 <- ggplot() +
  geom_point(data = air, 
             aes(x = Wind, y = Temp, 
                 color =  as.factor(Month)), 
             alpha = .6) +
  geom_line(data = prediction_viz, 
            aes(x = Wind, y = fit, 
                color = as.factor(Month))) +
  geom_ribbon(data = prediction_viz, 
              aes(x = Wind, y = fit, 
                  ymin = lwr, ymax = upr, 
                  group = as.factor(Month)), 
              alpha = .1) +
  scale_color_discrete(name = "Mois", 
                      labels = c("Mai", "Juin", "Juillet", 
                                 "Août", "Septembre")) +
  labs(title = "Température en fonction du vent",
       y = "Température prédite\n(avec intervalle de 95%)",
       x = "Vitesse du vent") +
  theme_minimal()

Visualiser des prédictions: Montrer les vraies observations

graph_pred2

Visualiser des prédictions: Montrer les vraies observations

1. J’ai sécifié une fonction ggplot vide;
2. J’ai tracé les points originaux, par couleur de mois (geom_point);
3. J’ai tracé les lignes de valeurs prédites, par couleur de mois (geom_line);
4. J’ai tracé les intervalles de confiance (geom_ribbon).

Faire la gigue: jitter

p = ggplot(data = air, aes(x = as.factor(Month), y = Temp)) +
  geom_boxplot() 

p + geom_point(alpha = .4)

Faire la gigue: jitter

Déplacer un peu des points, à l’horizontale et à la verticale.

p + geom_jitter(alpha = .4) # width = 0.2

Inclure des annotations: annotate

• Spécifier le type d’élément à inclure: du texte sur le graph, des lignes, des rectangles, des flèches, etc.
• Spécifier sa locatisation

p = ggplot(data = air, aes(x = Wind, y = Temp)) +
  geom_point() +
  theme_minimal()

p

Inclure des annotations: annotate

Des encadrés (rect)

p + annotate("rect", 
             xmin = 13, xmax = 16, 
             ymin = 87, ymax = 94,
             alpha = .2)

Inclure des annotations: annotate

Des traits (segment)

p + annotate("segment", 
             x = 19, xend = 15.5, 
             y = 95, yend = 91.5,
             color = "red",
             arrow = arrow(length = unit(.2,"cm")))

Inclure des annotations: annotate

Du texte (text)

p + annotate("text", 
             x = 17.5,
             y = 95,
             label = "blah blah blah", 
             size = 3)

Faire pivoter la banque de données

• Une banque de données peut être organisée dans un format “long” ou dans un format “large”
• Exemple au tableau

dat = data.frame(personne = c("A", "B", "C", "D"),
                 avant  = c(150, 115, 160, 133),
                 apres = c(152, 111, 161, 140)) 

Faire pivoter la banque de données: pivot_longer

Transformer au format “long”

# de la variable "avant" à la variable "apres"
dat_long = dat %>% 
  pivot_longer(avant:apres,
               names_to = "timing", 
               values_to = "mesure")

Faire pivoter la banque de données: pivot_longer

Transformer au format “long”

# toutes les variables sauf "personne"
dat_long = dat %>% 
  pivot_longer(!personne,
               names_to = "timing", 
               values_to = "mesure")

Faire pivoter la banque de données: pivot_wider

Transformer au format “large”

dat_large = dat_long %>% 
  pivot_wider(names_from = timing,
              values_from = mesure)