Exercice 02

On continu sur les données de baseball en testant cette fois ci le lien linéaire existant et en mettant en avant le fléau de la dimensionalité
Régression linéaire
Sélection automatique
Author

Clément Poupelin

Published

February 17, 2025

Modified

March 19, 2025

Intervenant.e.s

Rédaction

Relecture

Setup

Show the code 
# Données
library(ISLR)         # Hitters data 
library(dplyr)        # manipulation des données
Show the code 
set.seed(140400)

Données

On étudie toujours le jeu de données Hitters disponible dans la libraire {ISLR} de R. Il s’agit d’un jeu de données de la Major League Baseball provenant des saisons de 1986 et 1987.

Le jeu de données possède 322 lignes/individus pour les différents joueurs et 20 variables.
Parmi les variables, on trouve les informations suivantes :

AtBat Number of times at bat in 1986
Hits Number of hits in 1986
HmRun Number of home runs in 1986
Runs Number of runs in 1986
RBI Number of runs batted in in 1986
Walks Number of walks in 1986
Years Number of years in the major leagues
CAtBat Number of times at bat during his career
CHits Number of hits during his career
CHmRun Number of home runs during his career
CRuns Number of runs during his career
CRBI Number of runs batted in during his career
CWalks Number of walks during his career
League A factor with levels A and N indicating player's league at the end of 1986
Division A factor with levels E and W indicating player's division at the end of 1986
PutOuts Number of put outs in 1986
Assists Number of assists in 1986
Errors Number of errors in 1986
Salary 1987 annual salary on opening day in thousands of dollars
NewLeague A factor with levels A and N indicating player's league at the beginning of 1987

Comme pour l’Exercice 1, on va commencer par se débarasser des variables manquantes.

Show the code 
Hitters_Without_NA <- Hitters %>% na.omit()

Puis cette fois ci nous allons dans un premier temps nous concentrer sur un sous jeu de données composé des 18 premières lignes sans valeurs manquantes.

Show the code 
Hitters_Without_NA_18 <- Hitters_Without_NA[1:18, ]
Hitters_Without_NA_18 %>% dim()
[1] 18 20
Warning

Attention, on peut remarquer ici que le nombre de variable est supérieur au nombre d’individus. On est donc dans un cas classique de grandes dimension avec \(p>n\).

Maintenant, il conviendrait dans ce genre de situation d’effectuer de premières analyses descritptives. Mais celle ci ayant déjà été faite sur le jeu de données complet pendant l’Exercice 1, on se permettra de ne pas les refaires.

Analyse inférentielle

Modèle brut

On désire modéliser le salaire Salary en fonction des variables disponibles.

On va donc ajuster un modèle de régression linéaire en utilisant toutes les variables à disposition et analyser la qualité de cet ajustement.

Show the code 
mod1 <- lm(formula = Salary ~ .,
           Hitters_Without_NA_18) 
mod1 %>% summary()

Call:
lm(formula = Salary ~ ., data = Hitters_Without_NA_18)

Residuals:
ALL 18 residuals are 0: no residual degrees of freedom!

Coefficients: (2 not defined because of singularities)
             Estimate Std. Error t value Pr(>|t|)
(Intercept) -223.7909        NaN     NaN      NaN
AtBat         -3.2428        NaN     NaN      NaN
Hits          13.1990        NaN     NaN      NaN
HmRun        -60.8834        NaN     NaN      NaN
Runs           0.6875        NaN     NaN      NaN
RBI           10.3993        NaN     NaN      NaN
Walks          7.0114        NaN     NaN      NaN
Years         -2.3702        NaN     NaN      NaN
CAtBat         0.2643        NaN     NaN      NaN
CHits         -1.7919        NaN     NaN      NaN
CHmRun         5.3897        NaN     NaN      NaN
CRuns          4.0162        NaN     NaN      NaN
CRBI          -4.0134        NaN     NaN      NaN
CWalks         1.5822        NaN     NaN      NaN
LeagueN      233.6380        NaN     NaN      NaN
DivisionW    299.1771        NaN     NaN      NaN
PutOuts       -0.1250        NaN     NaN      NaN
Assists       -0.8539        NaN     NaN      NaN
Errors             NA         NA      NA       NA
NewLeagueN         NA         NA      NA       NA

Residual standard error: NaN on 0 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:    NaN 
F-statistic:   NaN on 17 and 0 DF,  p-value: NA

Résultats

On peut clairement constater que ce modèle brut ne fonctionne pas avec pourtant un \(R^2 = 1\). On retrouve donc le problème typique de l’analyse en grande dimension lorsque \(p>n\) (fléau de la dimensionalité).

On peut aussi s’amuser à regarder les critères AIC et BIC de ce modèles qui théoriquement se retrouve à tendre vers l’infini.

Show the code 
cat( "AIC = ", AIC(mod1), "et BIC = ", BIC(mod1))
AIC =  -Inf et BIC =  -Inf

Prediction

On va maintenant tenter de prédire la variable Salary pour les autres joueurs.
Déjà on peut regarder sur les 18 joueurs si la prédiction via le modèle nous donne des bonnes valeur.

Show the code 
Salary_hat <- predict(mod1, Hitters_Without_NA_18)
Salary <- Hitters_Without_NA_18$Salary
  • \(\widehat{Salary^{(1:18)}} - Salary^{(1:18)} =\) 0

Ce que l’on constate c’est qu’effectivement nous sommes avec un résultat qui pourrait nous faire penser que le modèle est bien ajusté avec une prédiction quasiment égale à la variable à prédire.

Pourtant si nous regardons la prédiction obtenue par le modèle pour les autres joueurs et que nous effectuons la même soustraction pour comparer la qualité de prediction, nous voyons bien l’inéfficacité du modèle.

Show the code 
Hitters_Without_NA_No18 <- Hitters_Without_NA[19:nrow(Hitters_Without_NA),]
Salary_hat_No18 <- predict(mod1, Hitters_Without_NA_No18)
Salary_No18 <- Hitters_Without_NA_No18$Salary
  • \(\widehat{Salary^{(\neg 1:18)}} - Salary^{(\neg 1:18)} =\) -70.88

En effet on voit bien au dessus que les valeurs ne sont en moyennes pas proche de 0.

Modèles parcimonieux

On va maintenant mettre un oeuvre une méthode de sélection automatique classique pour réduire le nombre de variable explicative et tenter d’éviter les problèmes de grande dimension.

Pour cela nous allons donc partir du plus petit modèle (celui avec seulement l’intercept) puis faire grandir le nombre de variable. Il va donc s’agir d’une méthode de sélection automatique forward.

Show the code 
mod0 <- lm(Salary~1, Hitters_Without_NA_18)
mod_forw <- step(mod0,
                 scope = formula(mod1),
                 trace = FALSE,
                 direction = c("forward"))
mod_forw %>% summary()

Call:
lm(formula = Salary ~ CWalks + League, data = Hitters_Without_NA_18)

Residuals:
    Min      1Q  Median      3Q     Max 
-215.51  -82.67  -48.10   26.13  302.49 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  259.3539    77.1270   3.363  0.00427 ** 
CWalks         0.9699     0.1606   6.039 2.27e-05 ***
LeagueN     -137.2850    79.1236  -1.735  0.10322    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 160.8 on 15 degrees of freedom
Multiple R-squared:  0.7495,    Adjusted R-squared:  0.7161 
F-statistic: 22.44 on 2 and 15 DF,  p-value: 3.095e-05

Résultats

Nous obtenons maintenant un modèle avec 2 variable dont une significative. Puis nous pouvons constater des valeurs assez élevés pour le \(R^2\) et \(R^2_{adjusted}\).

Et on a AIC = 238.692 et BIC = 242.253.

Donc sans aller tester si c’est un bon modèle prédictif, on constate déjà qu’il va s’agir d’un modèle descriptif fonctionnel avec \(n<p\)

Permutations

Maintenant, nous allons permuter de façon aléatoire les salaires des 18 joueurs et refaire la même analyse inférentielle. Ainsi, le lien linéaire devrait disparaitre et nous donner de mauvais résultats.

Note

pour des raisons de repouctibilité, une graine ou seed a été défini dans le setup afin que la génération aléatoire reste identique.

Faisons à nouveau le modèle brute sur nos 18 joueurs.

Show the code 
Hitters_Without_NA_18$Salary_permute <- sample(Salary)

mod1_permute <- lm(Salary_permute~., subset(Hitters_Without_NA_18, select = -Salary))
mod1_permute %>% summary()

Call:
lm(formula = Salary_permute ~ ., data = subset(Hitters_Without_NA_18, 
    select = -Salary))

Residuals:
ALL 18 residuals are 0: no residual degrees of freedom!

Coefficients: (2 not defined because of singularities)
              Estimate Std. Error t value Pr(>|t|)
(Intercept) -3424.4890        NaN     NaN      NaN
AtBat          38.2468        NaN     NaN      NaN
Hits          -95.2285        NaN     NaN      NaN
HmRun        -416.5695        NaN     NaN      NaN
Runs           40.5868        NaN     NaN      NaN
RBI            50.2860        NaN     NaN      NaN
Walks         -23.5877        NaN     NaN      NaN
Years         691.4211        NaN     NaN      NaN
CAtBat        -13.2054        NaN     NaN      NaN
CHits          53.0461        NaN     NaN      NaN
CHmRun         68.4706        NaN     NaN      NaN
CRuns         -16.2821        NaN     NaN      NaN
CRBI          -24.4509        NaN     NaN      NaN
CWalks          7.5470        NaN     NaN      NaN
LeagueN       304.4377        NaN     NaN      NaN
DivisionW    -445.0630        NaN     NaN      NaN
PutOuts         0.3775        NaN     NaN      NaN
Assists        -4.0206        NaN     NaN      NaN
Errors              NA         NA      NA       NA
NewLeagueN          NA         NA      NA       NA

Residual standard error: NaN on 0 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:    NaN 
F-statistic:   NaN on 17 and 0 DF,  p-value: NA

Résultats

A nouveau on peut constater l’inéfficacité d’un modèle avec toutes les variables du fait d’avoir \(p>n\).

Utilisons maintenant la sélection automatique en testant à nouveau l’approche forward.

Show the code 
mod0_permute <- lm(Salary_permute~1, subset(Hitters_Without_NA_18, select = -Salary))
mod_forw_permute <- step(mod0_permute, 
                         scope = formula(mod1_permute),
                         trace = FALSE,
                         direction = c("forward"))
mod_forw_permute %>% summary()

Call:
lm(formula = Salary_permute ~ Walks + HmRun, data = subset(Hitters_Without_NA_18, 
    select = -Salary))

Residuals:
    Min      1Q  Median      3Q     Max 
-437.54 -163.95    8.98  145.42  510.31 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  682.171    121.221   5.628 4.81e-05 ***
Walks        -13.212      5.056  -2.613   0.0196 *  
HmRun         22.524     15.187   1.483   0.1587    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 257.4 on 15 degrees of freedom
Multiple R-squared:  0.3583,    Adjusted R-squared:  0.2727 
F-statistic: 4.187 on 2 and 15 DF,  p-value: 0.0359

Résultats

On constate que plusieurs variables son significatives. Pourtant, on trouve ici un modèle avec de très mauvais \(R^2\) et \(R^2_{adjusted}\). Donc un modèle de mauvaise qualité avec en plus une variance assez grande.

Pour finir, on va maintenant reprendre le jeu de données Hitters complet et permuter tous les salaires de façon aléatoire. Ensuite, on va ajuster le meilleur modèle de régression possible pour expliquer les salaires en fonction des autres variables.

Show the code 
Hitters_Without_NA$Salary_permute <- sample(Hitters_Without_NA$Salary)

mod0_permute <- lm(Salary_permute~., subset(Hitters_Without_NA, select = -Salary))
mod1_permute <- lm(Salary_permute~1, subset(Hitters_Without_NA, select = -Salary))


mod_permute_back <- step(mod1_permute,
                         scope = formula(mod1_permute),
                         trace = FALSE,
                         direction = c("backward"))



mod_permute_forw <- step(mod0_permute,
                         scope = formula(mod1_permute),
                         trace = FALSE,
                         direction = c("forward"))


mod_permute_both <- step(mod0_permute,
                         scope = formula(mod1_permute),
                         trace = FALSE,
                         direction = c("both"))
Show the code 
mod_permute_back %>% summary()

Call:
lm(formula = Salary_permute ~ 1, data = subset(Hitters_Without_NA, 
    select = -Salary))

Residuals:
   Min     1Q Median     3Q    Max 
-468.4 -345.9 -110.9  214.1 1924.1 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   535.93      27.82   19.27   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 451.1 on 262 degrees of freedom
Show the code 
mod_permute_forw %>% summary()

Call:
lm(formula = Salary_permute ~ AtBat + Hits + HmRun + Runs + RBI + 
    Walks + Years + CAtBat + CHits + CHmRun + CRuns + CRBI + 
    CWalks + League + Division + PutOuts + Assists + Errors + 
    NewLeague, data = subset(Hitters_Without_NA, select = -Salary))

Residuals:
    Min      1Q  Median      3Q     Max 
-717.32 -311.48  -80.55  210.07 1729.22 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept)  213.83924  128.25879   1.667   0.0968 .
AtBat          0.79917    0.89573   0.892   0.3732  
Hits          -0.39046    3.35916  -0.116   0.9076  
HmRun        -15.65359    8.76187  -1.787   0.0753 .
Runs           1.54702    4.21144   0.367   0.7137  
RBI            3.30928    3.67472   0.901   0.3687  
Walks         -1.50102    2.58345  -0.581   0.5618  
Years         34.59958   17.53688   1.973   0.0496 *
CAtBat        -0.19099    0.19107  -1.000   0.3185  
CHits          0.76040    0.95306   0.798   0.4257  
CHmRun         2.69347    2.28496   1.179   0.2396  
CRuns         -0.99130    1.06030  -0.935   0.3508  
CRBI          -1.04738    0.97858  -1.070   0.2855  
CWalks         0.99760    0.46354   2.152   0.0324 *
LeagueN     -188.25478  111.98651  -1.681   0.0940 .
DivisionW     59.97471   57.03348   1.052   0.2940  
PutOuts       -0.05287    0.10941  -0.483   0.6294  
Assists       -0.33618    0.31253  -1.076   0.2831  
Errors        -0.51571    6.20483  -0.083   0.9338  
NewLeagueN   240.36156  111.62089   2.153   0.0323 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 445.9 on 243 degrees of freedom
Multiple R-squared:  0.09395,   Adjusted R-squared:  0.0231 
F-statistic: 1.326 on 19 and 243 DF,  p-value: 0.1672
Show the code 
mod_permute_both %>% summary()

Call:
lm(formula = Salary_permute ~ AtBat + HmRun + Years + CAtBat + 
    CWalks + League + Assists + NewLeague, data = subset(Hitters_Without_NA, 
    select = -Salary))

Residuals:
   Min     1Q Median     3Q    Max 
-681.1 -325.6 -101.4  220.0 1725.7 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  238.07953  116.47262   2.044  0.04198 *  
AtBat          0.90098    0.28464   3.165  0.00174 ** 
HmRun         -7.42307    4.29005  -1.730  0.08479 .  
Years         35.62737   15.88548   2.243  0.02578 *  
CAtBat        -0.14624    0.04353  -3.360  0.00090 ***
CWalks         0.64126    0.25014   2.564  0.01094 *  
LeagueN     -195.98139  109.10111  -1.796  0.07363 .  
Assists       -0.33430    0.22402  -1.492  0.13686    
NewLeagueN   245.27466  108.53924   2.260  0.02468 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 439.3 on 254 degrees of freedom
Multiple R-squared:  0.0807,    Adjusted R-squared:  0.05174 
F-statistic: 2.787 on 8 and 254 DF,  p-value: 0.005669

Résultats

On constate qu’aucune méthode de sélection de variable ne permet d’avoir ne serait-ce qu’un modèle correct ce qui montre bien qu’avec la permutation aléatoire de la variable Salary, le lien linéaire qui existait à disparu.

Conclusion

Dans un premier temps, on a pu avoir un aperçu de ce qu’il se passe lorsque l’on se retrouve face au fléa de la dimensionalité avec un sous jeu de données où le nombre de variables était supérieur au nombre d’individus.

Puis, on a aussi pu voir l’importance du lien linéaire dans la construction d’un modèle de régression. Cela renforce par l’exemple la véracité du modèle de régression linéaire (au cas où l’on en doutais encore).

Session info

Show the code 
sessioninfo::session_info(pkgs = "attached")
─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.4.2 (2024-10-31)
 os       Ubuntu 24.04.1 LTS
 system   x86_64, linux-gnu
 ui       X11
 language (EN)
 collate  fr_FR.UTF-8
 ctype    fr_FR.UTF-8
 tz       Europe/Paris
 date     2025-03-19
 pandoc   3.2 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/x86_64/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
 package    * version date (UTC) lib source
 dplyr      * 1.1.4   2023-11-17 [1] CRAN (R 4.4.2)
 forcats    * 1.0.0   2023-01-29 [1] CRAN (R 4.4.2)
 ggplot2    * 3.5.1   2024-04-23 [1] CRAN (R 4.4.2)
 ISLR       * 1.4     2021-09-15 [1] CRAN (R 4.4.2)
 kableExtra * 1.4.0   2024-01-24 [1] CRAN (R 4.4.2)
 lubridate  * 1.9.4   2024-12-08 [1] CRAN (R 4.4.2)
 purrr      * 1.0.4   2025-02-05 [1] CRAN (R 4.4.2)
 readr      * 2.1.5   2024-01-10 [1] CRAN (R 4.4.2)
 stringr    * 1.5.1   2023-11-14 [2] CRAN (R 4.3.3)
 tibble     * 3.2.1   2023-03-20 [2] CRAN (R 4.3.3)
 tidyr      * 1.3.1   2024-01-24 [1] CRAN (R 4.4.2)
 tidyverse  * 2.0.0   2023-02-22 [1] CRAN (R 4.4.2)

 [1] /home/clement/R/x86_64-pc-linux-gnu-library/4.4
 [2] /usr/local/lib/R/site-library
 [3] /usr/lib/R/site-library
 [4] /usr/lib/R/library

──────────────────────────────────────────────────────────────────────────────